{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2a0b97d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp2d, interpn\n",
    "from numpy.matlib import repmat\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f6be88",
   "metadata": {},
   "source": [
    "# Entry/Exit with Homogeneous firms\n",
    "## and translog demand\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "244c912c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sequence_jacobian import solved, simple, create_model\n",
    "\n",
    "@simple\n",
    "def household(w, psi, nu):\n",
    "    L = (w/psi)**nu\n",
    "    return L\n",
    "\n",
    "@solved(unknowns={'v': 1.0, 'N': 1.0}, targets=['valfunc', 'accumulation'])\n",
    "def firms(C, sigma, v, beta, delta, Ne, N, L):\n",
    "    mu      = 1 + 1/(sigma*N)\n",
    "#    y       = p**(-sigma)*C\n",
    "    ell     = L/N            # imposing labor market clearing\n",
    "    pi      = (1-1/mu)*(C/N)\n",
    "    valfunc = v - (pi + beta*(1-delta)*v(1))\n",
    "    accumulation =  N - (1-delta)*N(-1) - Ne\n",
    "    return mu, ell, pi, valfunc, accumulation\n",
    "    \n",
    "@simple \n",
    "def mkt_clearing(mu, N, v, fE, w, Z, C, L, pi, Ne, Nbar):\n",
    "    p       = mu*w/Z\n",
    "    rho     = math.e**(-0.5*(Nbar - N)/(2*Nbar*N))\n",
    "    goods_mkt = p - rho                   # as in BGM\n",
    "    free_entry = v - fE*w/Z\n",
    "    budget     = C - w*L - N*pi\n",
    "    return goods_mkt, free_entry, budget, rho\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a70a25",
   "metadata": {},
   "source": [
    "### Steady State"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d3fbc144",
   "metadata": {},
   "outputs": [],
   "source": [
    "ces = create_model([household, mkt_clearing, firms], name=\"CES Model\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "bb885ca8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4190476190476196"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "calibration = {'beta':    0.96,\n",
    "               'psi':     1,\n",
    "               'fE':      .25,\n",
    "                'nu': 2.0,\n",
    "                'delta':  0.11,\n",
    "                'sigma':  .35,\n",
    "                'w': 1,\n",
    "                'Z': 1.01,\n",
    "                'Ne': .75, \n",
    "                'C': .25, \n",
    "                'Nbar': 1}\n",
    "\n",
    "ss = ces.steady_state(calibration)\n",
    "ss['mu']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "283c6fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "w0    = .15\n",
    "Nbar0 = 1/(calibration['sigma']*(1/w0-1))\n",
    "Ne0   = calibration['delta']*Nbar0\n",
    "\n",
    "calibration['Ne'] = Ne0\n",
    "\n",
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'Nbar': Nbar0, 'w': w0, 'Ne': Ne0}, \n",
    "      targets={'rho': 1.0, 'free_entry': 0.0, 'goods_mkt': 0.0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1c603fd4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.37644679813266"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss['budget']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "87cc0f0d",
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "No convergence after 100 iterations",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-20-5c837566a03e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m ss=ces.solve_steady_state(calibration, \n\u001b[0m\u001b[1;32m      2\u001b[0m       \u001b[0munknowns\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0;34m'w'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'w'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Nbar'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Nbar'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Ne'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mss\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Ne'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'C'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;36m.05\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m       targets={'free_entry': 0.0, 'goods_mkt': 0.0, 'budget': 0.0, 'rho': 1.0})\n",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/blocks/block.py\u001b[0m in \u001b[0;36msolve_steady_state\u001b[0;34m(self, calibration, unknowns, targets, dissolve, options, **kwargs)\u001b[0m\n\u001b[1;32m    157\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mcompute_target_values\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtargets\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mss\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    158\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 159\u001b[0;31m         _ = solve_for_unknowns(residual, unknowns, solver, opts['solver_kwargs'],\n\u001b[0m\u001b[1;32m    160\u001b[0m                                \u001b[0mtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mopts\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'ttol'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mopts\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'verbose'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    161\u001b[0m                                \u001b[0mconstrained_method\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mopts\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'constrained_method'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/blocks/support/steady_state.py\u001b[0m in \u001b[0;36msolve_for_unknowns\u001b[0;34m(residual, unknowns, solver, solver_kwargs, residual_kwargs, constrained_method, constrained_kwargs, tol, verbose)\u001b[0m\n\u001b[1;32m    201\u001b[0m         \u001b[0;31m# If no bounds were provided\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mbounds\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 203\u001b[0;31m             unknown_solutions, _ = solvers.broyden_solver(residual_f, initial_values,\n\u001b[0m\u001b[1;32m    204\u001b[0m                                                           tol=tol, verbose=verbose, **solver_kwargs)\n\u001b[1;32m    205\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/utilities/solvers.py\u001b[0m in \u001b[0;36mbroyden_solver\u001b[0;34m(f, x0, y0, tol, maxcount, backtrack_c, verbose)\u001b[0m\n\u001b[1;32m    117\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Too many backtracks, maybe bad initial guess?'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    118\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 119\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'No convergence after {maxcount} iterations'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    120\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    121\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: No convergence after 100 iterations"
     ]
    }
   ],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'w':2*ss['w'], 'Nbar': ss['Nbar'], 'Ne': ss['Ne'], 'C': .05}, \n",
    "      targets={'free_entry': 0.0, 'goods_mkt': 0.0, 'budget': 0.0, 'rho': 1.0})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52073d0e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.471767566296194e-13"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss['mu']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e52b213",
   "metadata": {},
   "source": [
    "# Shock to Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d08066a",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 300\n",
    "dZ = -.03 * .685 ** np.arange(T)\n",
    "irf = {}\n",
    "irf = ces.solve_impulse_linear(ss, ['w', 'Ne', 'C'], ['free_entry', 'goods_mkt', 'budget'], {'Z': dZ})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c00e2675",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABjHklEQVR4nO3deXhU1fnA8e8bliRASFgDhH0X2WVz3xWtilr3Fa0/6lrtYlurttba2mpr61qr1Yy7orKILAKKKFqQHcISlrAkgSSEkJB9fX9/3Js4hgQGmGSWvJ/nmSdz5t65887kzJ33nnvuOaKqGGOMMcaEk4hAB2CMMcYY42+W4BhjjDEm7FiCY4wxxpiwYwmOMcYYY8KOJTjGGGOMCTuW4BhjjDEm7FiC04BE5C4RyRSRAhHpEKAYeouIikjzepY/JiJvu/d7urE2O8I2J4vIkoaI1xhfHKleG2NMSCY4IrJTRMpEpGOtx9e4O73eAQrNO5YWwDPABaraRlX3+2m7DZZcqOpuN9bKhth+Q3D/3/0DHYepXyh8X41pbO7BZPWtSkSKvco3ugef5bXW+7X73C9FpMR9LFtEpolI10C/p2ATkgmOawdwfXVBRIYB0YEL5xDxQBSwIdCBGBME/PZ9tVYbEw7cg8k2qtoG2A1c6vXYO+5qH3ivp6pPeW3iXve5A4E44J+N+w6CXygnOG8Bt3iVbwXe9F5BRH4kIqtF5KCIpIrIY17LokTkbRHZLyK5IrJcROLdZZNFJEVE8kVkh4jcWFcAIhIpIv8SkT3u7V/uYwOBZHe1XBH5op7nTxCRb93XXysiZ3ktOyQGETkBeBk42c3cc4/0Pr3c7sa4V0R+WU88P2j2P9LnICJ/F5ED7rKLvB7/UkSecN9bgYjMEpEOIvKOG+Ny76N2ERksIgtEJEdEkkXkGq9lHhF5UURmu3EsE5F+7rKv3NXWuq9zrYh0FJFP3c80R0S+FpFQrufh4rDf1yN8V6vr5U9EZDdwyPdJRH7sthQNdevME17LzhKRNK/yThF5SEQ2uvU3UUSi/Px+jWkUqpoDfAwMDXQswSaUd/xLgbYicoI4fUauBd6utU4hzk41DvgRcJeIXO4uuxWIBXoAHYA7gWIRaQ08B1ykqjHAKcCaemJ4GJgAjARGAOOAR1R1C3Ciu06cqp5T+4kikgDMBp4A2gO/Aj4WkU71xaCqm9w4/+dm83E+vM9qZwMDgAuA34rIefW8p+r4jvQ5jMdJ4joCTwGviYh4Lb8OuBlIAPoB/wMS3fe6CfiD1+ssAN4FOuMc5b8kIid6bet64I9AO2Ab8GcAVT3DXT7C/Tw+AH4JpAGdcFrRfgfYfCSBd6Tvqy91+EzgBOBC7wdF5Dbgb8B5qprkYzw3utvph3ME/MjRvBljgoU4p35/DKwOdCzBJpQTHPj+qPB8YDOQ7r1QVb9U1fWqWqWq64D3cHaSAOU4iU1/Va1U1ZWqetBdVgUMFZFoVd2rqvWdZroReFxVs1R1H86P8M0+xn4TMEdV57jxLQBWABcfZQxHep/V/qiqhaq6HifRuP6QDR3qcDHsUtVX3f46bwBdcRKKaomqul1V84C5wHZVXaiqFcCHwCh3vUuAnaqaqKoVqroK52jkKq9tTVPV79znvoOTUNan3I2ll6qWq+rXahOuBYt6v68+1uHH3Dpc7PXYA8CDwFmquu0oYnlBVVPdo98/49v3wZjGdo3bGl196+a17Dm3FX8tsBf4RUAiDGLhkODcAEym1ukpABEZLyKLRGSfiOThtH509HruZ8D77qmbp0SkhaoW4hxd3gnsdU+NDK7n9bsBu7zKu9zHfNELuNq78gKnAV2PMoYjvc9qqUcTpw8xZHitW+TebeO1PNPrfnEd5ep1ewHja30ONwJd6notoKjW69T2NE4rz3z39NpvD7OuaVz1fl+PoQ5XexB4UVXT6lh2OEf1fTAmQKaqapzXbY/Xsp+5jyWo6o3uQbbxEtIJjqruwum8eDEwrY5V3gU+AXqoaixO/xVxn1uuqn9U1SE4p18uwe0joKqfqer5OC0Bm4FX6wlhD84PdLWe7mO+SAXeqlV5W6vqX48QQ12tEfW+Ty89jjbOo/gcjkcqsLjW59BGVe86lo2par6q/lJV+wKXAr8QkXP9GrE5Jkf4vvpSh+uq+xcAj4jIj70eKwRaeZW7cKij/j4YY0JLSCc4rp8A57gtDrXFADmqWiIi43COHgEQkbNFZJjbH+AgzqmNShGJF5HL3L4hpUABUN9l0+/h7Fw7uedBf8+h/YDq8zZwqYhcKCLNxOn0fJaIdD9CDJlAdxFp6cv79PKoiLRy+7bcBnxwuOCO8nM4Hp8CA0XkZhFp4d7GitOh2heZQN/qgohcIiL93f5AB92YQ+ay9yagvu+rL3W4LhuAicCLInKZ+9ga4GIRaS8iXXBOY9V2j/tda4/TT+uw3wdjTOgJ+QTH7eexop7FdwOPi0g+TvIx1WtZF+AjnB/BTcBinKQjAqej6h4gB6cfwN31bP8JnH4z64D1wCr3MV/iTgUm4exc9+G0ZDzovv7hYvgCZ6eeISLZPrzPaotxTt18DvxdVecfIcSj+RyOmarm4xyFX+e+VgZOh9FIHzfxGPCGe3rrGpyO1AtxErL/AS+p6pd+Dtsco8N8X32pw/Vtcy1OC+yr4lzN9xZOv4SdwHzqTl7edZeluDefvrfGmNAh1v/SGNOUiMhO4A5VXRjoWIwxDSfkW3CMMcYYY2o7YoIjIqe6/TAQkZtE5BkR6XWk5xljjGl6RGSiOAN2brOrGE0gHfEUlYiswxnEbjjOue3XgCtVtfYYFcYYY5ow96KNLThjHaUBy4HrVXVjQAMzTZIvp6gq3IHSJgHPquqzOFc8GGOMMd7GAdtUNUVVy4D3cX47jGl0vkxaly8iD+GMvHuGm6G3aNiw6taxY0ft3bt3IF7ahImVK1dmq2qnQMdxJFbXzfEKUF1P4IeDKKbhTOtSr5iYGO3RowetWjlDF+3fv59WrVoRHR2NqpKTk3NIuXXr1kRFRVFVVcWBAwcOKbdp04bIyEgqKyvJzc09pBwTE0PLli2pqKggLy/vkHLbtm1p0aIF5eXlHDx48JBybGwszZs3p6ysjPz8/EPKcXFxNGvWjNLSUgoKCg4pt2vXjoiICEpKSigsLDyk3L59e0SE4uJiioqKDil36NABgKKiIoqLi+stFxYWUlpaSvv27essFxQUUF5eTrt27eos5+fnU1lZSVxcXJ3lgwcPoqrExsbWWwZo27YtAHl5eYhIveXc3FyaNWtGTExMneUDBw7QokUL2rRp84NycnJy3XVdVQ97w7mc+hfA6W65J3DLkZ7XELeTTjpJjTkewAr1Y50EXgeygKR6lgvOnF7bcIYTGO3Ldq2um+Pl77ruyw24GvivV/lm4Pk61puCM8TGioSEBF29erWqqlZUVGhiYqKuXbtWVVXLyso0MTFR169fr6qqxcXFmpiYqBs3blRV1cLCQk1MTNTNmzerqmp+fr4mJibq1q1bVVU1NzdXExMTdfv27aqqmpOTo4mJibpjxw5VVd23b58mJibq7t27VVU1MzNTExMTNS0tTVVV9+7dq4mJibp3715VVU1LS9PExETNzMxUVdXdu3drYmKi7tu3T1VVd+zYoYmJiZqTk6Oqqtu3b9fExETNzc1VVdWtW7dqYmKi5ufnq6rq5s2bNTExUQsLC1VVdePGjZqYmKjFxcWqqrp+/XpNTEzUsrIyVVVdu3atJiYmakVFhaqqrl69WhMTE2v+5ytWrNA33nijpvzdd9/p22+/XVP+3//+p++++25N+ZtvvtEPPvigpvz111/rhx9+WFP+8ssv9eOPP64pf/HFFzpjxoya8oIFC/STTz6pKX/22Wf66aef1pTnzp2rc+fOrSl/+umn+tlnn9WUP/nkE12wYEFNecaMGfrFF1/UlD/++GP98ssva8offvihfv311zXlDz74QL/55pt667ovfXBaAyWqWinOLNmDgbmqWn7YJ/pARCYCzwLN3C/FXw+3/pgxY3TFivqGvDHmyERkpaqO8eP2zsAZc+dNVT1kNl8RuRi4D2f03vE4p3kPe0QLVtfN8fN3XffxNU/GmTPsQrf8EICqPlnfc6yum+NVX133pQ/OV0CkOLNff44zCq7HDwE1A14ELgKGANeLyJDj3a4xjUlVv8IZCLE+k3CSH1XVpUCciHQ91tdLzy0+8krGBM5yYICI9HFHW78OZwoOY45ZZZVysOTo21R8SXBEnckUr8RparwCOPGoX+lQR90Zbf/+/axZswaAyspKPB4P69atA6C8vByPx0NSUhIAJSUleDweNm3aBDjnJj0eD8nJyYBzrtHj8bBtmzMBcV5eHh6Ph5SUFMA5t+fxeNi5cycA2dnZeDweUlOd08tZWVl4PB7S050JkTMyMvB4PGRkOPNCpqen4/F4yMrKAiA1NRWPx0N2tjP48M6dO/F4PBw4cACAlJQUPB4PeXl5AGzbtg2Px0NBQQEAycnJeDweioqceS03bdqEx+OhpKQEgKSkJDweD+XlTiVYt24dHo+HykpnloI1a9bg8XhqPsuVK1fy5pvfz3e4fPly3nnnnZry0qVLee+992rK3377LVOnfj+47JIlS/joo49qyosXL2batO+nF1q0aBEzZ86sKS9cuJBZs2bVlOfPn8/s2bNryvPmzWPevHk15dmzZzN//veDLc+aNYuFC78fl23mzJksWrSopjxt2jQWL15cU/7oo49YsmRJTXnq1Kl8++23BEBdfRIS6lpRRKaIyAoRWbFv36Hz5n22IYMznlrEip2Hy6eMCRxVrQDuxZnIeBPOZJEbAhuVCWaqSl5xOZv2HuSLzZm8tXQXf5u3mQfeX801L/+P0/72BYMemcv/vXH0rXy+dDIWt9nxRpx5ZMA5pXS8fOqMJiJTcM7XkpBQ5++CMcGs9oSRUPekkajqK8Ar4DTb115+Wv+OdI6J5JEZSXx632k0b2bjdJrgo6pzgDmBjsMEh8oqJf1AMWm5RezJLWFvbjF78opJr76fW0xh2Q+nC2weIXSNi6JrbDRje7ena2wUg7oc/cXbvvTBOQP4FfCNqv5NRPoCD6jqz4761X643auBC1X1Drd8MzBOVe+r7zl2rtYcr4bolyAivYFP6+mD8x/gS1V9zy0nA2ep6t7DbbO+uj53/V7uemcVv79kCLef1scv8ZvwFIg+OMfC9uvhoapKSc8tZktmPlsyC9y/+WzLKqC0ouoH63ZsE0m3uCi6xUbTNS6KhLhousZG082937FNJBERdR0b1q2+un7EFhy3j8FXXuUU4LiSG1ca0MOr3B1nskVjwsknwL0i8j5OC2XekZKbw5k4tAtnDOzEMwu2cMnwrnRuG+W3QI0x5khUlb15JSRn5rPVTWa2ZuazNauAIq+WmK6xUQyIj+Hkvh0YEN+GHu1b0S02mi6xUUS18MdJoCPz5RRVQ6npjAak43RGuyGA8ZgQUVWlZBeWsie3hD1uE2e6+3dPbgktm0fw8V2nNEosIvIecBbQUUTSgD/gjhOlqi/jNNVfjHOZeBFOJ/3jeT3+eNmJXPjPr3hi9iaeu37U8WzOGGPqVVZRxbq0XNam5bE1M5/kzHy2ZRaQX1pRs06nmEgGxrfhmjE9GNQlhoHxbejfOYbY6IAMl/cDAUtwVLVCRKo7ozUDXrfOaAagsLSCve452roSmIy8Esoqf9jk2bplMxLaRdMtLpo+HVs3Wqyqev0Rlitwjz9fs0/H1tx5Zl+e+2Ib143twSn9O/pz88aYJqq0opK1qXksS9nP0h37WbnrACXlzr62XasWDIyP4fJRCQzsEsPAzm0YGB9Du9YtAxx1/QLZgmOd0Zq43KIykjPy2ZJVwJYM9+ggq4CcwrIfrBch0KVtFN3iohnZI45uw6JJiHPK1be2Uc0R8f2cbai7++z+TF+TzqMzk5h7/xm0bG4djo0xR6ekvJK1qbksTclhacp+Vu0+UNNfZnCXGK4b25MJfdszulc7OrWJDLl97BETHPcU0n1Ab+/1VfWyhgvLhJP8knK2ukmMd+ezrPzSmnViIpszsEsMF54YT8/2rZ0OaG7yEh8TaVcM1RLVohl/vOxEbves4LUlO7jrrH6BDskYE+RKyitZvTuXpSn7WbZjP6t251JWUYUInNClLTeO78X4vu0Z17t9ULfM+MqXFpwZODOIzwKqDr+qacrKKqrYkpnvtsrk1yQ03oPTRbdoxoD4NpwxsBMD450mzkFdYujSNirkjg4C7ZzB8Zw/JJ7nPt/KZSO7kRAXHeiQjDFBpLJKWbEzh2+272dpyn7WpDoJTYTAkG5tuWVCL8b37cC43u2JbRX4PjP+5kuCU6KqzzV4JCbklFVUsTYt1zlfm5LDil05NedrWzaLoG+n1ozp3Y4b4ns6iUx8DN3bRR/V5X/m8P5w6RDOe2Yxf5q1kZdvPinQ4RhjgsCWzHymrUpnxup0Mg6WECEwNCGWW0/uxYS+HRjTu31QdAJuaL4kOM+KyB+A+UDNOQVVXdVgUZmgdLgOaNXna8f0bsfgLm3p3aGVnVZqBN3bteK+cwbw9GfJfJmcxVmDOgc6JGNMAOzLL+WTtXuYtiqNDXsO0ixCOGtgJx655ATOGNiJtlHhn9DU5kuCMwxnRthz+P4UlbplE8ZKKypZs9vpgLbMTWhK3fO1g7u05fpxPZngNm+Gw/naUHXH6X34eGUaf/hkA5890KHRxpgwxgRWSXkl8zdmMn1VGl9tzaayShmWEMsfLh3CpSO60bFNZKBDDChfEpwrgL7ufFEmjJWUV7Im1e2AlpJT06PeuwPahL7tGdenPXGtLKEJFpHNm/H4pKHc9Noy/rM4hfvPGxDokIwxDaSqSlm2I4fpq9OYsz6DgtIKusVG8dMz+nLl6AT6dz76KQ3ClS8JzlogDshq2FBMIFRWKd9uz2b66nTmJWVQVFaJCJzYrS03TehV00ITjh3QwslpAzryo+FdefHLbVw+qhu9OjTeWEDGmIa3LauA6avTmLF6D+m5xbRu2YyLhnXlytEJTOjTwfo21sGXBCce2Cwiy/lhHxy7TDyEbc44yPRV6cxYk07mwVJiIptz6fBunDcknnF9mkYHtHDz6I+G8OXmLB77ZAOvTx5rV6UZE+IOlpQzfVU601alsTYtjwiB0wd04tcTB3HBkC5Et7TT0YfjS4LzhwaPwjSKrIMlzFyzh2mr09m09yDNI4QzB3bi0UsSOO+EeOu7EeK6xEbx8/MH8sTsTczfmMmFJ3YJdEjGmGNQUl7J20t38cKibeQWlTOka1se+dEJXDaim80/dxR8mWxzsYjEA2Pdh75TVTtdFSKKyir4bEMG01al8822bKoURnSP5TG3E1qHJt4JLdzcekpvPlyRxuOzNnL6gI60ahnQwcqNMUehskqZsTqdZxZsIT23mDMGduJXFwxkePe4QIcWknwZyfga4GngS0CA50XkQVX9qIFjM8eopl/NqnTmbXD61STERXP3Wf25YnQC/Tq1CXSIpoG0aBbBny4fyjX/+R8vfLGNX08cHOiQjDFHoKp8mbyPv83bzOaMfIZ3j+Xpq4bbPHPHyZfDu4eBsdWtNiLSCVgIWIITZA7pVxPVnMtGdOOKUQmM7d3eOqE1EeP6tOfK0Qm8+nUKV47uTv/OltAaE6xW7z7AX+duZtmOHHp3aMWLN4zm4mFdrA+dH/iS4ETUOiW1H7AR3ILIurRcnpqXzJJt2TSPEM4a1InfX9Kdc0/obP1qmqiHLjqBBRsz+cMnSbz9k/G2szQmyGzfV8DT85KZtyGDjm0i+dPlQ7lubA9a2ACpfuNLgjNPRD4D3nPL12IzgAeF7fsK+Mf8ZOasz6BdqxY8dNFgrjqpu/WrMXSKieTBCwfx+5kb+HTdXi4d0S3QIRljgMyDJfxr4VamrkglqnkEvzh/ID85rQ+tI62/nL8d9hMV57DvOZwOxqfh9MF5RVWnN0Jsph5784p5duFWPlyZRlTzCO4/dwB3nN6HmCY4FLep343jezF1RSpPzN7I2YM708Z2oMYEzMGScv6zeDuvLdlBZZVy84Re3HtO/yY/2nBDOuweT1VVRGao6knAtEaKydTjQGEZL325jTf+twsUbjm5F/ecbV8QU7dmEcITlw/jipe+4V8LtvDIJUMCHZIxTU5pRSVv/e/7S74njezGL88fRM8OrQIdWtjz5ZBuqYiMVdXlDR6NqVNhaQWvL9nBK1+lUFhWwRWjuvPAeQPo0d6+IObwRvaI47qxPUn8didXjenO4C5tAx2SMU2CqvLJ2j08NS+Z9NxiTh/Qkd9MHMzQhNhAh9Zk+JLgnA38VER2AYU4p6lUVYc3aGSGsooq3vtuN89/sZXsgjIuGBLPry4cxMB4m2vE+O7XFw5iXtJeHp2RxNSfnmwdjo1pYKUVlfxh5gbeX57K0IS2/O3HwzltgF3y3djqTXBEpI+q7gAu8veLisjVwGPACcA4VV3h79cIZZVVysw1zmBPaQeKGd+nPa/cMpjRPdsFOjQTgtq1bslvLxrMbz5ez8er0rnqpO6BDsmYsJWRV8Kdb69kTWou95zdj1+cP4hmNkRHQByuBecj4CTgdVU918+vmwRcCfzHz9sNaarK55uyePqzZJIz8zmxW1v+fMUwzhjQ0Y66zXG5+qQefLA8lSfnbOL8E+Jt8lRjGsCKnTnc9c4qCksr+PeNo7loWNdAh9SkHS7BiRCRPwADReQXtReq6jPH+qKqugmwH20vK3fl8Jc5m1m56wB9OrbmhRtGcfHQrjY4n/GLiAjhT5cP5dLnl/DHTzfwj6tH2PfPGD9RVd5Ztps/ztpAQlw079wx3roSBIHDJTjXAZe76wTsPyUiU4ApAD179gxUGA2mqkr59+Lt/GN+Mp1iIvnLFcO4ekx3G+zJ+N2J3WK575wBPPv5Vnq1b8395w0IdEjGhDzv/jZnDerEs9eOshbSIFFvgqOqycDfRGSdqs492g2LyEKgrumMH1bVmb5uR1VfAV4BGDNmjB5tHMEsr7icX05dw8JNWVw6oht/vXKYDfZkGtQD5w0g7UAx/1y4ha5xUVwzpkegQzImZHn3t7n37P78/PyB1t8miPgym/hRJzfu8847luc1FRv25HHX26vYk1vMY5cO4dZTetspA9PgRIS//ngYWfklPDRtPZ1jIjlrUOdAh2VMyFm+M4e73l5FUVkFL980molDrb9NsLHzIAHw4YpUrnzpW8oqqvjgpxOYfGofS25Mo2nRLIJ/33QSg+JjuPudVSSl5wU6JGNChqry1tJdXP/KUtpENmPGPadachOkApLgiMgVIpIGnAzMdue6Cnsl5ZU8NG0dD360jpN6tePTn53GSb3aBzos0wS1iWyO57axtGvVksmJy0nNKQp0SMYEvdKKSn778XoenZHE6QM6MvPe06wzcRDzKcERkVNE5AYRuaX6djwvqqrTVbW7qkaqaryqXng82wsFqTlFXPXyt7z3XSp3n9WPt34y3qZYMAHVuW0Ub9w+lvLKKm5N/I4DhWWBDsmEMBG5WkQ2iEiViIwJdDz+lpFXwrX/WcoHK1K59+z+/PfWscRGW2fiYHbEBEdE3gL+jjPZ5lj3FnaVtyEtSs7ikueXsGt/Ea/eMoZfTxxsHdFMUOjfOYb/3jqGtAPF3PHmCkrKKwMdkgld1eObfRXoQPxt+c4cLnl+CVsz83n5ptH86kIbvC8U+HLJzhhgiKqG1RVMjaGySnn28608/8VWBndpy8s3jaZXh9aBDsuYHxjbuz3/unYk97y7igfeX8OLN462nbc5auE4vpmq8vbSXfxx1kZ6tG/Fe/83ngF2Sipk+HKKKom6L/c2h5FTWMbkxO947vOtXDmqO9PuOsWSGxO0Lh7WlUd/NIR5GzL406cbseMZ05BEZIqIrBCRFfv27Qt0OHUqrajkNx+v49GZGzh9QEdm3HOqJTchxpcWnI7ARhH5DiitflBVL2uwqELc2tRc7n5nFfvyS3nyymFcN7ZHWB3VmPB0+2l92JNbzH+X7KBbXBRTzugX6JBMkGkq45tVVSm/+GAts9fv5b5z+vPz8wbaqPIhyJcE57GGDiJcVA/X/fisjXSKieSju05mePe4QIdljM9+d/EJ7D1Ywl/mbKZLbDSXjegW6JBMEGkq45v96/OtzF6/l99eNJg7z7REP1T5MtDfYhGJx+lcDPCdqmY1bFihp7iskoenr2fa6nTOHNiJf107knatWwY6LGOOSkSE8I+rR7Avv5RfTV1LpzaRnNyvQ6DDMqbRzFyTznOfb+Xqk7rz0zP6Bjoccxx8uYrqGuA74GrgGmCZiFzV0IGFkr15xVzx0jdMX5POz88bSOLksZbcNCEiMlFEkkVkm4j8to7lZ4lInoiscW+/D0Scvopq0YxXbx5Drw6tmPLWCpIz8gMdkgkB4TC+2erdB3jwo3WM692eJ64Yal0LQpwvnYwfBsaq6q2qegswDni0YcMKHXnF5Ux+fTlpB4pJnDyW+88bYOdqmxARaQa8CFwEDAGuF5Ehdaz6taqOdG+PN2qQxyC2VQs8t4+jVctmTE78jr15xYEOyQS5UB/fbE9uMVPeWkl820j+fdNoIps3C3RI5jj5kuBE1Doltd/H54W90opKfvrWClKyC3jl5pNsTp+maRywTVVTVLUMeB+YFOCY/CIhLprEyePIL6ngtsTlHCwpD3RIxjSIorIK7nhjBcVllbx261g62CCsYcGXRGWeiHwmIpNFZDIwG5jTsGEFv6oq5cEP17E0JYenrxrBKf07BjokExgJQKpXOc19rLaTRWStiMwVkRMbJ7TjN6RbW16+6SS2ZRVw51srKauoCnRIxvhVVZXywPtr2JxxkOdvGGVTL4SRIyY4qvogzuV8w4ERwCuq+puGDizY/e2zzXyydg+/mTiYy0fV9Xtmmoi6zkfWvux1FdBLVUcAzwMz6txQkI4NctqAjjx11XC+3b6fX3+0lqqqoLuq15hj9vf5yczfmMkjPxrC2dYKH1Z8uUwcVf0Y+LiBYwkZb3y7k/8sTuHmCb2480zrZd/EpQE9vMrdgT3eK6jqQa/7c0TkJRHpqKrZtdYL2rFBrhzdnb15JTz9WTLxsVH8duJg64BpQt60VWm89OV2rh/Xk9tO7R3ocIyf1ZvgiMgSVT1NRPL54RGpAKqqbRs8uiA0LymDx2Zt4Pwh8Tx22Ym2kzfLgQEi0gdIB64DbvBeQUS6AJmqqiIyDqfldH+jR3qc7j6rH3vzivnP4hSy88t44vKhRLe0jpgmNK3clcNvP17PhL7teXyS7cvDUb0Jjqqe5v61E5KulbtyuP/91YzsEcdz142y+XoMqlohIvcCnwHNgNdVdYOI3Okufxm4CrhLRCqAYuC6UJzbTUT442VD6dA6kue+2MrGvQf5942j6d3RpiAxoSXtQBFT3lxJt7goXr7pJFo0s+tmwpGvs4kf8bFwt31fAT95YwXd4qJ57daxduRqaqjqHFUdqKr9VPXP7mMvu8kNqvqCqp6oqiNUdYKqfhvYiI9dswjh5+cP5PXJY9mTW8ylLyxh/oaMQIdljM8KSiv4iWcFZZVVvDZ5LHGtbMyycOVL2vqDKz5EpDlwUsOEE5yy8ku49fXvaB4hvHHbONrbIH6miTt7UGc+ve80endozZS3VvK3eZupqLQrrExwq6xS7n9vNdv2FfDSjaPp16lNoEMyDajeBEdEHnL73wwXkYPuLR/IBHyeVC3UFZRWcLtnOfsLynjt1rH07NAq0CEZExR6tG/Fh3eezPXjevLvL7dzy+vfkV1QeuQnGhMgf5u3mc83Z/HYpUM4fUCnQIdjGli9CY6qPun2v3laVdu6txhV7aCqDzVijAFTXlnFPe+sYtPefF66cTQjesQFOiRjgkpUi2Y8eeUwnr5qOCt3HeBHz33Nyl05gQ7LmENMXZ7KK1+lcMvJvbj55N6BDsc0Al/GwXlIRNqJyDgROaP6djwvKiJPi8hmEVknItNFJO54ttcQVJXfTVvP4i37+PPlQzl7sI2PYEx9rh7Tg2l3n0Jk82Zc+5+lJH6zgxDsR23C1LKU/Tw8Yz2nD+jI7y+payYVE4586WR8B/AVzlUif3T/Pnacr7sAGKqqw4EtQNC1CP1r4VY+XJnGz84dwHXjegY6HGOC3ondYpl132mcNagTf5y1kfveW01haUWgwzJN3K79hdz59kp6tG/FCzeMprldMdVk+PKfvh8YC+xS1bOBUcBxDbOqqvNVtXrPtxRncLSg8f53u3n2861cfVJ3fn7egECHY0zIiI1uwSs3j+HXEwcxZ/1eJr34DduybDZyExgHS8r5yRsrqFJ4/daxxEa3CHRIphH5kuCUqGoJgIhEqupmYJAfY7gdmOvH7R2XRZuzeHhGEmcM7MRfrhxmgz8Zc5QiIoS7z+rPWz8Zz4HCMia98A2frttz5Cca40cVlVXc++5qdmYX8vJNJ9l4TU2QLwlOmttHZgawQERmUmso+rqIyEIRSarjNslrnYeBCuCdw2yn0ebnWZeWy93vrOKErjG8dONoG/zJmONwav+OfPqz0xjUJYZ7313N47M2Um6XkptG8uTczXy1ZR9PXD6Uk/t1CHQ4JgCOOBeVql7h3n1MRBYBscA8H5533uGWi8itwCXAuYcb1bWx5ufZvb+I2z3L6dCmJa9PHkubSJ+m6TLGHEbX2Gjen3Iyf5mzide/2cG6tFxeuGE0XWKjAh2aCWNJ6Xm8/s0Obp7Qy/pQNmGHGwenrfu3ffUNWA8sAY5rdCQRmQj8BrhMVYuOZ1v+kFNYxq2J31FRpXhuG0fnGNv5GuMvLZtH8NhlJ/LsdSPZsOcglzz/NZ9tyLCrrEyDUFX+PHsTcdEteHCiP3tTmFBzuHMw77p/VwIr6vh7PF4AYnBOea0RkZePc3vHTFW5771VpOcW899bxtC/s41saUxDmDQygZn3nkq7Vi356Vsr+fG/v2VZSsjNOWqC3KLkLP6Xsp/7zx1A2yjrVNyUHW6yzUvcv338/aKq2t/f2zxWM9fs4Ztt+3ni8qGM6d0+0OEYE9YGxscw5/7T+WhlGv9auIVrX1nKWYM68eCFgzixW2ygwzMhrqKyiifnbKZPx9bcML5XoMMxAebLODgzReR6EQm7OQryist5YvYmRvSI4wY7T2tMo2jRLILrx/Vk8YNn89BFg1m9O5cfPbeEn73nXPFizLGauiKNrVkF/GbiYFo2t4tEmjpfasAzwOnAJhH5UESuEpGw6KTyzwVb2F9YyhOThhIRYZeDG9OYolo046dn9uOrX5/NPWf3Y8HGTM57ZjEPT19P5sGSQIdnQkxhaQXPLNjCmF7tuPDE+ECHY4KAL1M1LFbVu4G+OFczXQNkNXRgDS0pPY83/7eTm8b3Ylh3axo3JlBio1vw4IWDWfzrs7hhfE8+WJ7KmU8v4m/zNpNXVB7o8EyI+M9XKWQXlPLwj06w8csM4FsLDiISDfwYuBNnVOM3GjKohlZVpTw6M4l2rVryqwusl70xwaBzTBSPTxrK5788k4knduHlxds5/akv+PeX2ykuqwx0eCaIZR4s4dWvUrhkeFdG9WwX6HBMkPClD84HwCbgHOBFoJ+q3tfQgTWkD1emsnp3Lr+7+ARiW1kve2OCSa8OrfnXdaOYfd/pjOndnr/N28yZTy/i7aW7bKBAU6dn5m+hoqqKX184ONChmCDiSwtOIk5Sc6eqfqGqIb2HOVBYxl/nbmZc7/ZcOToh0OEYY+oxpFtbXp88lg/vPJme7VvxyIwkzntmMZ+s3UNVlY2hYxyb9h5k6spUbj25Nz07hN21MOY4+JLgfAU8JCKvAIjIABG5pGHDajhPfZbMwZIKHr/8RDtPa0wIGNu7PR/eeTKvTx5DdItm/Oy91Vz07Nf89+sUsvKtM3JT9+TczbSNasG95wTN6CMmSPjaglMGnOKW04AnGiyiBrR69wHeX76b207pzeAubQMdjjHGRyLCOYPjmfOz0/nXtSNp2TyCJ2Zv4uQnv2By4nfMXJNu/XSaoK+27OOrLfu475z+xLVqGehwTJDxZcKlfqp6rYhcD6CqxRKCTR+VbsfizjGRPHD+wECHY4w5BhERwuWjErh8VALbsvKZtiqdGavTuf/9NbSJbM7EoV24clQCE/p2sKEfwlxllfKXOZvo0T6am0+2Qf3MoXxJcMrcq6gUQET6AaUNGlUDeGfZLpLSD/L89aNsIk1jwkD/zjH8euJgfnXBIJbtyGH66jTmrs/go5VpdI2NYtLIBK4cncDA+JhAh2oawMer0tickc/z148isnmzQIdjgpAvv/R/wJk9vIeIvAOcCkxuyKD8bV9+KU9/lsxp/TtyyfCugQ7HGONHERHCyf06cHK/Djw+aSgLNmYyfXU6r36dwsuLt3Nit7ZcMSqBy0Z2s4l0w0RxWSX/mJ/MyB5xtk839TpigqOqC0RkFTABEOB+Vc1u8Mj86Mm5mygpr+SPk6xjsTHhLKpFMy4d0Y1LR3Qju6CUWWv3MH11Ok/M3sRf5mzi9AGduHJ0AhcM6UJ0Szvq9ycReRq4FKfP5nbgNlXNbYjX+u/XKWQeLOWFG0bbPt3Uq94ER0RG13por/u3p4j0VNVVDReW/3y3I4dpq9K55+x+9OtkM4Ub01R0bBPJbaf24bZT+7Atq4Dpq9OYsXoP97+/htYtm3FK/45M6NuB8X3ac0LXtjSzPjvHawHwkKpWiMjfgIeA3/j7Rfbll/Ly4u1ceGI8Y22CZHMYh2vB+Yf7NwoYA6zFacEZDiwDTmvY0I5feWUVj85IIiEumnvPHhDocIwxAdK/cxsevHAwvzx/EN/tzGHmmj18uz2bBRszAWgb1ZxxfTowoW97JvTtYAnPMVDV+V7FpcBVDfE6/1q4hdKKKn4z0Qb1M4dXb4KjqmcDiMj7wBRVXe+WhwK/apzwjs8b3+4kOTOfV24+yZqjjTFERAgT+nZgQt8OAOzNK2ZZSg5LU/azNGU/Czc5CU9MVHPG9W5fs+6QbpbwHKXbgQ/qWygiU4ApAD179vR5o9uy8nl/eSo3T+hFX2uRN0fgSyfjwdXJDYCqJonIyIYLyT8y8kr454ItnDO4M+cPsZlljTGH6hobXXPZOTj7jWU7nGRnWUoOn2925hWOiWzO2D7ta1p4hnRtS/NmPk3lF1ZEZCHQpY5FD6vqTHedh4EK4J36tqOqr+BM3syYMWN8Hpb6r3M306pFM352rrXImyPzJcHZJCL/Bd7GuVT8Jpy5qYLan2ZvpKJKeexS61hsjPFNF/fy8kkjnYQn82CJ27qTw7Id+/nCTXjaRDZnVM84BsXHMLBLDIPiYxgQ34ZWLcN7CApVPe9wy0XkVuAS4FxV9et8Gt9uz2bhpix+M3Ew7VvboH7myHz5Nt4G3AXc75a/Av7dYBH5wZKt2cxet5efnzfQ5iYxxhyz+LY/THiyDpawdEcOy1L2syY1l7eW7qK04vvp+Xq0j3aTHSfpGRgfQ99OrYlqEf6nyEVkIk6n4jNVtcif265yB/VLiIvmtlN7+3PTJoz5cpl4CfBP9xb0Sisq+f3MJHp3aMVPz+wb6HCMMWGkc9soLhvRjctGdAOc0XR35xSRnJHP1sx8kjPz2ZpZwJfJ+6hwJwSNEOjdsfUPEp9BXdrQq0NrWoTXaa4XgEhggdtqvlRV7/THhj9Zu4ek9IP889oRTSJZNP4RkPZUEfkTMAmoArKAyaq6xx/b/u/XO0jJLuSN28fZF8EY06CaRQh9OramT8fWTBz6fdeUsooqdu4vZEtmPlsynMQnOSOfzzZkUD0ReotmQtfYaLrFRdEtLpqEuGi6ubeEuCi6xkbTOoRGXVfVBpntsqS8kqc/S2ZoQlsmjUhoiJcwYSpQ356nVfVRABH5GfB74Lgz/dScIp7/YisXDe3CmQM7He/mjDHmmLRsHsFA9xQVw79/vKS8km1ZBWzNymdLZgHpB4rZk1vM0u37yThYUpP8VIuNblGT8HTzSoC6xTrlzjGRYd/ZOfGbnaTnFvP01cNtfjFzVI6Y4IjI1ar64ZEeOxqqetCr2Bp3nqvj9finG4kQ4dFLhvhjc8YY41dRLZoxNCGWoQmxhyyrqKwiM7+UPbnF7q2k5n7agWK+25HDwZKKHzynb8fWfPGrsxop+saXU1jGS4u2ce7gzpzSr2OgwzEhxpcWnIeA2slMXY8dFRH5M3ALkAecfZj1fBov4fNNmSzYmMlvLxpMt7jo4wnNGGMaXfNmESS4p6rqk19Szt68EtLdxKdZmF8h+tznWykqr+Shi21QP3P0DjdVw0XAxUCCiDzntagtzhgHh3Wk8RJU9WHgYRF5CLgXZ1LPQ/gyXkJJeSWPzdpA/85tuP3UPkcKzRi/cq8eeRZoBvxXVf9aa7m4yy8GinD6nIXEVCcmuMREtSAmqkWTmCE9ZV8Bby/dxbVje9C/c/i/X+N/h2vB2QOsAC4DVno9ng/8/EgbPtJ4CV7eBWZTT4Lji5cWbSM1p5j3/m8CLZuH9/loE1xEpBnwInA+kAYsF5FPVHWj12oXAQPc23icYRbGN3asxoSSp+YlE9k8ggfOs0H9zLGpNxtQ1bWq+gYwDHhbVd9wyzOB0uN5URHxrrGXAZuPdVs7sgt5eXEKk0Z24+R+HY4nLGOOxThgm6qmqGoZ8D7OFYLeJgFvqmMpECciXQ+30f3797NmzRoAKisr8Xg8rFu3DoDy8nI8Hg9JSUkAlJSU4PF42LTJGX+zqKgIj8dDcnIyAAUFBXg8HrZt2wZAXl4eHo+HlJQUAA4cOIDH42Hnzp0AZGdn4/F4SE1NBSArKwuPx0N6ejoAGRkZeDweMjIyAEhPT8fj8ZCV5QyCl5qaisfjITs7G4CdO3fi8Xg4cOAAACkpKXg8HvLy8gDYtm0bHo+HgoICAJKTk/F4PBQVOUOpbNq0CY/HQ0lJCQBJSUl4PB7Ky8sBWLduHR6Ph8rKSgDWrFmDx+Op+SxXrlzJm2++WVNevnw577zz/SC7S5cu5b333qspf/vtt0ydOrWmvGTJEj766KOa8uLFi5k2bVpNedGiRcycObOmvHDhQmbNmlVTnj9/PrNnz64pz5s3j3nz5tWUZ8+ezfz530/jNGvWLBYuXFhTnjlzJosWLaopT5s2jcWLF9eUP/roI5YsWVJTnjp1Kt9++y2hbPnOHOZtyODOM/vROSYq0OGYEOVLc8d8wPukcDSwsJ51ffVXEUkSkXXABXw/iOAxOX1ARx6++ITjDMmYY5IApHqV09zHjnYdRGSKiKwQkRXVP97GNEWRzSM474TO3HG6jWVmjp0caTRtEVmjqiOP9FhjGDNmjK5YsaKxX9aEERFZqapj/Li9q4ELVfUOt3wzME5V7/NaZzbwpKouccufA79W1ZV1bROsrpvj5++63lCsrpvjVV9d96UFp1BERntt6CSg2J/BGRPC0oAeXuXuOP3XjnYdY4wxfuTLZeIPAB+KSPUOuStwbYNFZExoWQ4MEJE+QDpwHXBDrXU+Ae4VkfdxOhfnqerexg3TGGOaliOeogIQkRbAIECAzaoakA4CIrIP2FXP4o5AdiOGcywsRv84nhh7qapfh7kWkYuBf+FcJv66qv5ZRO4EUNWX3cvEXwAm4lwmfpuqHrZN/jB1PRT+PxAacYZ7jH6v6w3B6nqjCPcY66zrvvTBiQLuBk7DGXH4a+BldxLOoCEiK4L9fLPF6B+hEGNDCZX3HgpxWozBLVTeeyjE2VRj9OUU1Zs4Y98875avB94CrvZnIMYYY4wx/uJLgjNIVUd4lReJyNqGCsgYY4wx5nj5chXVahGZUF0QkfHANw0X0jF7JdAB+MBi9I9QiLGhhMp7D4U4LcbgFirvPRTibJIx+tIHZxNOB+Pd7kM9gU1AFaCqOtzfQRljjDHGHA9fEpxeh1uuqvVd1WSMMcYYExA+XSZujDHGGBNKQmrqbRGZKCLJIrJNRH5bx3IRkefc5eu8R2BupPh6iMgiEdkkIhtE5JA5tkTkLBHJE5E17u33jRmjVxw7RWS9G8MhY7IEwWc5yOszWiMiB0XkgVrrBMVn2RCsrvs1VqvrQczqut/iDOp67sbQuHVdVUPihjOI2nagL9ASWAsMqbXOxcBcnAEJJwDLGjnGrsBo934MsKWOGM8CPg2Cz3Mn0PEwywP6Wdbxv8/AGcwp6D7LBnq/Vtf9F6vV9SC9WV33a5whU8+9/vcNWtdDqQVnHLBNVVNUtQx4H5hUa51JwJvqWArEiUjXxgpQVfeq6ir3fj5OZ+xDZo0OEQH9LGs5F9iuTae/l9X1xmV1PXCsrjeeYKrn0Ah1PZQSnAQg1aucxqGVzJd1GoWI9AZGAcvqWHyyiKwVkbkicmLjRlZDgfkislJEptSxPGg+S5z5nd6rZ1kwfJb+ZnXdv6yuBy+r6/4TSvUcGqGu+zLQX7CQOh6r3UPal3UanIi0AT4GHlDVg7UWr8JpkisQZw6jGcCARg4R4FRV3SMinYEFIrJZVb/yWh4sn2VL4DLgoToWB8tn6W9W1/3L6nrwsrruPyFRz6Hx6nooteCkAT28yt2BPcewToMSZ2LSj4F3VHVa7eWqelBVC9z7c4AWItKxMWN0X3uP+zcLmI7TVOwt4J+l6yJglapm1l4QLJ9lA7C67kdW14Oa1XU/CaF6Do1U10MpwVkODBCRPm72dx3wSa11PgFucXuLTwDyVHVvYwUoIgK8BmxS1WfqWaeLux4iMg7nf7C/sWJ0X7e1iMRU3wcuAJJqrRbQz9LL9dTTjBkMn2UDsbruJ1bXg57Vdf/EGEr1HBqprofMKSpVrRCRe4HPcHpfv66qG0TkTnf5y8AcnJ7i24Ai4LZGDvNU4GZgvYiscR/7Hc7oz9UxXgXcJSIVQDFwnao2djNhPDDdrUPNgXdVdV6QfZaISCvgfOCnXo95xxgMn6XfWV33K6vrQczqut+ERD2Hxq3rNtCfMcYYY8JOKJ2iMsYYY4zxiSU4xhhjjAk7luAYY4wxJuxYgmOMMcaYsGMJjjHGGGPCjiU4xhhjjAk7luAYY4wxJuxYgmOMMcaYsGMJjjHGGGPCjiU4xhhjjAk7luAYY4wxJuxYgmOMMcaYsGMJTpATERWR/oGOw5imQER+JyL/DXQcxpjjZwnOMRKRnSJSLCIFXrcXAh1XoInIYyLydqDjMP4jIjeIyAq3ju8Vkbkiclqg4zpeInKWiKR5P6aqf1HVOwIVkwk/IvKQiMyp9djWeh67rnGjC2+W4ByfS1W1jdft3kAHZIw/icgvgH8BfwHigZ7AS8CkAIZlTCj5CjhVRJoBiEgXoAUwutZj/d11jZ9YguNnIjJZRL4RkX+KSK6IpIjIKe7jqSKSJSK3eq3vEZGXRWSBiOSLyGIR6VXPtmNF5E0R2Sciu0TkERGJEJFIEckRkWFe63Z2W5g6VR+pisiv3dffKyKXi8jFIrLFfe7vvJ4bISK/FZHtIrJfRKaKSHt3WW/3tNmtIrJbRLJF5GF32UTgd8C17tH+2ob6nE3DE5FY4HHgHlWdpqqFqlquqrNU9UG33v1LRPa4t3+JSKT73Oo690uvOneb17YvFpGNbp1PF5FfuY9PFpElteKoOU3rfl9ecluRCtzvWhf3tQ+IyGYRGeX13J3uEfRGd3miiESJSGtgLtDNqwW2W+0WSBG5TEQ2uN/lL0XkhFrb/pWIrBORPBH5QESiGua/YULYcpyEZqRbPgNYBCTXemw7cKGIbHK/Fyki8lPvDbn78L3u9+2OWt+NSBH5u7tfznR/V6Ib4f0FLUtwGsZ4YB3QAXgXeB8Yi5Oh3wS8ICJtvNa/EfgT0BFYA7xTz3afB2KBvsCZwC3Abapa6r7GTV7rXg8sVNV9brkLEAUkAL8HXnXXPwk4Hfi9iPR11/0ZcLn7Gt2AA8CLtWI5DRgEnOs+9wRVnYdzpP+B26I14nAfkgl6J+PUmen1LH8YmICzkx4BjAMe8VreBae+JgA/AV4UkXbusteAn6pqDDAU+OIo4rrGfZ2OQCnwP2CVW/4IeKbW+jcCFwL9gIHAI6paCFwE7PFqgd3j/SQRGQi8BzwAdALmALNEpGWtWCYCfYDhwOSjeB+mCVDVMmAZThKD+/drYEmtx74CsoBLgLbAbcA/RWQ01BxA/gI4D+e35MxaL/U3nPo90l1eva9vulTVbsdwA3YCBUCu1+3/cHZwW73WGwYoEO/12H5gpHvfA7zvtawNUAn0cMuKU1mb4ezMh3it+1PgS/f+eCAViHDLK4Br3PtnAcVAM7cc4253vNe2VgKXu/c3Aed6LesKlAPNgd7uc7t7Lf8OuM69/xjwdqD/P3bzSx2/Ecg4zPLtwMVe5QuBne796jrX3Gt5FjDBvb/brb9ta21zMrCk1mMK9Hfve4BXvZbdB2zyKg8Dcr3KO4E7vcoXA9u9Ykyr9Vo19Rd4FJjqtSwCSAfO8tr2TV7LnwJeDvT/zW7Bd3Pr1XT3/lpgAE5i7P3YrXU8bwZwv3v/deBJr2X9+f73QYBCoJ/X8pOBHYF+74G8WQvO8blcVeO8bq+6j2d6rVMMoKq1H/NuwUmtvqOqBUAOTsuJt45AS2CX12O7cLJ0VHUZTgU/U0QG41T6T7zW3a+qld4x1RFndUy9gOlus3wuTsJTidMHo1qG1/2iWu/HhIf9QEcRaV7P8m4cWh+96+1+Va3wKnvXkx/jJBu73NOyJx9FXLXr7eG+W+D1/aojxsP5wftT1Sp3Wwle69j3wPjiK+A0twWzk6puBb4FTnEfGwp8JSIXichSt9tALs53pKO7jW78sC573+8EtAJWeu2357mPN1mW4ASHHtV33FNX7YE9tdbJxmlF8e6f0xPniLLaGzinnW4GPlLVkmOMJxW4qFbyFqWq6Ud8pnNEYcLD/4ASnNOVddnDofWxdr2tk6ouV9VJQGeco9Sp7qJCnB01UNP58nj18LrvHeOR6uoP3p+IiLstX74Hxnj7H87p2inANwCqehCnjk1x/+4BPgb+jtPiH4dzWlTcbewFuntt07teZ+Mk9yd67bNjVbVJJ9yW4ASHi0XkNPfc/p+AZarqnZ3jtr5MBf4sIjHidET+BeB9SfZbwBU4Sc6bxxHPy+7r9AIQp6Oyr1fNZAK9RcTqVohT1Tycc/gvup3SW4lIC/co8ymc/imPuPWjo7vuEYcIEJGWInKjiMSqajlwEKeFEJym+hNFZKTbYfcxP7yVe0Skuzgd5X8HfOA+ngl0cDtT12Uq8CMROVdEWgC/xDlN/K0fYjJNiKoW43Qb+AVO/5tqS9zHvsJpoY8E9gEVInIRcIHXulOB20TkBBFphVf/Grd18VWcPjudAUQkQUQubLh3FfzsR+j4zJIfjoNTX2fMI3kX+APOqamTcPo+1OU+nCPcFJwvxrs452UBUNU0nM6Wyg+/REfrWZzTW/NFJB9YitPHxxcfun/3i8iq44jBBAFVfQZnB/wIzo43FbgXp9XlCZyd9jpgPU7de8LHTd8M7BSRg8CduB3kVXULzpVbC4GtOPX8eL0LzMf53qRUx6iqm3GStBS3Wf8Hp65UNdmN63mcI+RLcYaGKPNDTKbpWYzTYuldp792H/tKVfNxLvCYinNhxw14dTNQ1bnAczhXYG3DaRUCJ+kG+I37+FL3e7UQ50KQJkvczkgmQETEg9PR8ZEjrevj9l7HuTLEL9szJpSJyE7gDlVdGOhYjPEnd8iCJCCyVl8347IWnDAiIr2BK3EuwTXGGBNGROQK9xRvO5zLwmdZclM/S3DChIj8CSebf1pVdwQ6HmOMMX73U5xTxdtx+q3dFdhwgpudojLGGGNM2LEWHGOMMcaEnfoG8ApKHTt21N69ewc6DBPCVq5cma2qQT/4ldV1c7ysrpumor66HtAEx51b41mcaQj+q6p/Pdz6vXv3ZsWKFY0SmwlPIrLryGsd9TYPW4/dAeKexRmVtAiYrKqHvYTe6ro5Xg1R1318Xduvm0ZVX10P2CkqcaaJfxFnwrshwPUiMiRQ8RhzLHysxxfhzD0zAGfU0n83apDGNBLbr5tgEsg+OOOAbaqa4g6c9T5w2NFy9+/fz5o1awCorKzE4/Gwbt06AMrLy/F4PCQlJQFQUlKCx+Nh06ZNABQVFeHxeEhOTgagoKAAj8fDtm3bAMjLy8Pj8ZCSkgLAgQMH8Hg87Ny5E4Ds7Gw8Hg+pqc4Aw1lZWXg8HtLTnVHbMzIy8Hg8ZGQ4U9Okp6fj8XjIysoCIDU1FY/HQ3Z2NgA7d+7E4/Fw4MABAFJSUvB4POTl5QGwbds2PB4PBQUFACQnJ+PxeCgqKgJg06ZNeDweSkqc2RiSkpLweDyUl5cDsG7dOjweD5WVzgCxa9aswePx1HyWK1eu5M03vx/sePny5bzzzveTmC9dupT33nuvpvztt98yderUmvKSJUv46KOPasqLFy9m2rRpNeVFixYxc+bMmvLChQuZNWtWTXn+/PnMnj27pjxv3jzmzZtXU549ezbz58+vKc+aNYuFC78fymTmzJksWrSopjxt2jQWL15cU/7oo49YsuT78bSmTp3Kt982yAC0vtTjScCb6lgKxIlI18Nt1Oq61fVqQVTXfWH7davrNeVA1/VAJjgJ/HCysDR+OIkdACIyRURWiMiK6n+yMUHEl3psdd00FVbXTdAI2GXiInI1cKGq3uGWbwbGqep99T1nzJgxaudqzfEQkZWqOsaP2ztiPRaR2cCTqrrELX8O/FpVV9a3Xavr5nj5u677+Jq2XzeNrr66HsgWnDR+OBtqd3ycidiYIOJLPba6bpoKq+smaAQywVkODBCRPu4s2tfhNbGYMUdDVckpLCMpPY+lKfsb86V9qcefALeIYwKQp6p7GzNIEz5UleyCUtal5bJiZ06gw6nN9uvGb6qqlKz8Etam5rJ694Gjfn7ALhNX1QoRuRf4DOdywtdVdUOg4jHBraS8kr15JezJLSY9t5g9NbcS9uQ590vKqwDo0jaKpb87t1Hiqq8ei8id7vKXgTk4l4hvw7lM/LZGCc6EpOKyypo67dT3EvbmFruPOd+B0gqnrg/o3IYFvzgzwBF/z/br5mgUllawN8+p43tyi9mb+/39PXnF7M0toazSqeuje8Yx7e5Tj2r7AR0HR1Xn4Oz8TRNXUVnFzv1FbMvKJ+2AdxLjVPb9hWWHPKdzTCTd4qI5oUtbzh3cmW5x0XSLiyYhLrpRY6+rHruJTfV9Be5p1KBM0CqvrGJndiFbswpIO1DEntwS0nOL2esmMDm16roIxMdE0S0uihO7teX8IfF0i42iW1w0Pdq3CtC7qJ/t1021sooqUrIL2JZVQNqBQxOY3KIfdjCPEIhv69Tt4d3jmDg0ioS4aLrFRtOrw9HX9ZAaydiEvqoqJe1AMVsy80nOzGdLZj5bMgvYnlVQk6kDtGrZjIS4aLrGRTM0oS3dYqN/kMDEx0YS2bxZAN+JMYdXWaWk5hQ59Twjny1ZBWzJyCclu4Dyyu8v7mjdshkJ7aJrduoJcdF0i4uqqfPxbaNo2dxm1THBq6Kyil05RU49zyyo2b/vzC6kour7uh4T1dyt39GM7hVXsz+v3rfHx0TSvJn/6rolOKZBqCoZB0tIzshna2YByZn5bHWTmeLyypr1usVGMbBLDKcP6MjA+BgGxrehV/vWtI1ujjMAsDHBTVVJzy2uqedb3NvWzIKaU0kA3dtFMzA+hrMHd2ZQlzYM6BxDj/ataBtldd2Ehqoqp64nZ+SzJctJ3JMzC9i+r4Ayt66LQI92rRgYH8MFQ+IZ1CWG/p3buHW9RaPGawmOOW6qSmpOMUtT9rM6NbdmB59fUlGzTqeYSAbGt+G6cT0YFB/DgPgYBsS3afQKb8zxUFV2ZBeybEcOa1Nz3cS9gILS7+t6fNtIBsbHcNOEXgyKj2Ggu4NvE2m7WxM6VJXt+wr4X0oO69z9+tasAorKfniAOiD+hweo/Tu3oVXL4KjrwRGFCSmqyu6cIpam7GdpSg7LUvazJ88ZeTM2ugWDusQwaWQ3Z+fu3tq1bhngqI05eqpKSnYhS1P2sywlh6Up+8nKLwWgXSunrv94dAID4mMY1CWGgZ1jiG1lSbsJParK1qwCllXv13fsJ7vA6Q/WsU1LBnWJ4ZoxPZx6HiIHqJbgmCNSVXbtr05o9rNsRw573YSmY5uWjO/bgbv6tGdC3w7079zGmttNyHKOWgtr6vnSlP3scxOazjGRTOjbwb21p0/H1lbXTciqqnITmh37axL46os5usZGccaATkzo24HxfdvTs32rkKzrluCYQ6gqO70SmqUp+8k86OzkO7aJZELf9ozv24GT+7anXydLaEzo8m6Gr97JZxc4dT2+bSSn9OtQk9T07hCaO3ljwElotmTls3S7k7wv25FTc8VeQlw0Zw5yEpoJfTrQo310WNR1S3AMAHnF5Xy2IYMlW7N/0AzfqeaotT3j+3SgXyc7ajWh7UBhGXOTMvhmW/YPmuG7tI3itP7fJzS9LKExIS67oJQ56/fy7bb9LNuxnwPuZdkJcdGcPagzE/o6Le/BONyAP1iC04SVVVSxeMs+pq9OY+GmLMoqqn7QDD++b3v6WjO8CQOlFZV8sSmLaavT+TI5i/JKrWmGH+/u5EO1Gd4Yb8VllSzYlMn0VWl8tTWbyiqle7tozj0h3tmv92kftglNbZbgNDGqytq0PKavSuOTtXs4UFRO+9YtuWFcT64YlcDw7rG2kzdhQVVZsesA01alM3vdHg6WVNApJpLJp/Tm8lEJDOna1uq6CQtVVcrSHfuZviqduUkZFJRW0C02iiln9OXKUU4n+KbIEpwmIjWniBmr05m+Op2U7EJaNo/g/CHxXDkqgTMGdqKFHwdXMiaQdmQXMn1VGtPXpJOaU0x0i2ZceGI8V4zuzqn9Ovh1IDFjAmlrZj7TVqczc3U6e/JKaBPZnIuGduGK0QlM6NOBiIimncBbghPG8orKmb1+L9NXp7F8pzNR2fg+7fnpmX25aFjXoL/Ezxhf5RSW8em6PUxblc6a1FxE4NR+HXng3IFMHNqF1jYGjQkT+/JL+WTtHqavTiMp/SDNIoQzBnTktxefwPknxBPd0kZ4r2bf+jDzg341G7Moq6yiX6fWPHjhICaN7Eb3dk3j3KsJfyXllXyxOYtpq5x+NRVVyuAuMTx00WAmjUygS2xUoEM0xi+KyyqZvzGD6avT+drtVzMsIZbfXzKES0d0o1NMZKBDDEqW4ISJHdmFJH6zg1luv5oOrVtyw/ieXDk6gWEJ1q/GhI8tmfkkfrOD2ev2crCkgs4xkdx2am+uGNWdId3aBjo8Y/wmKT0Pz7c7mWf9ao6JJTghLiOvhGc/38rUFak0ixDrV2PCVmpOEf9cuIXpq9OJat6MiUO7cMWoBE7t35FmTbyvgQkvKfsK+MeCLcxet9f61RwHS3BCVG5RGf9evB3PNzupUuXmCb245+z+1lRpwk52QSkvfLGNd5btQkT4v9P7cteZ/Wz6DxN2vA9YWzaL4L5z+vN/Z/S1/pLHyBKcEFNUVkHiNzt5efF2CkoruGJkAj8/f2CTGdfANB35JeW8+vUO/vt1CiXllVwzpgf3nzeArrHRgQ7NGL+qfcB60/ie3HvOADtgPU6W4ISIsooqPli+m2c/30Z2QSnnndCZX104iMFdrM+BCS8l5ZW8vXQXLy7axoGicn40rCu/uGAg/Tq1CXRoxviVHbA2LEtwglxVlTJr3R7+MX8Lu3OKGNe7Pf+5eTQn9Wof6NCaPBF5GrgUKAO2A7epam4d6+0E8oFKoEJVxzRimCGjorKKaavS+efCLezNK+H0AR158MJBDO8eF+jQjPErO2BtHAFJcHz9YWjKVJVFyVk8NS+ZzRn5DOnalsTbxnLWwE52RVTwWAA8pKoVIvI34CHgN/Wse7aqZjdeaKFDVZmXlMHf5yezfV8hI3rE8Y+rR3BK/46BDs0cBRG5GngMOAEYp6orAhtR8KmqUj5Zu4dnFtgBa2MIVAvO0fwwNDnLd+bw1LzNLN95gF4dWvHc9aO4ZFhX6z0fZFR1vldxKXBVoGIJVd9sy+apeZtZm5ZH/85tePmmk7jwxHhL4kNTEnAl8J9ABxJs7IA1MAKS4NgPQ9027T3I3z9L5vPNWXSKieSJy4dy7dgedrl3aLgd+KCeZQrMFxEF/qOqr9S1kohMAaYA9OzZs0GCDBZrU3N5+rNklmzLpltsFE9fNZwrR3e3y71DmKpuAuwHuxY7YA2cYOiDc7gfhiahqKyCP36ykakrU4mJbM5vJg5m8im9bcjtICAiC4EudSx6WFVnuus8DFQA79SzmVNVdY+IdAYWiMhmVf2q9kpu4vMKwJgxY9QvbyDI5BWX8+iMJD5Zu4f2rVvy6CVDuHF8T6JaWF034SWnsIzffryO+Rsz7YA1QBoswfHTD0PYH9Vu31fAXW+vZFtWAVNO78vdZ/UntpWNeRAsVPW8wy0XkVuBS4BzVbXOpERV97h/s0RkOjAOOCTBCXeb9h7kzrdXkn6gmJ+543vE2PgeIcWX/bqP2wnr/frq3Qe4551VZBeW8eCFg7j91D52wBoADZbg+OOHwd1O2B7Vzl2/lwc/WkfL5hG8eft4ThtgnSpDiYhMxOk7dqaqFtWzTmsgQlXz3fsXAI83YphB4eOVaTw8Yz2x0S14f8oExvS2TpWh6Ej79aPYTlju11WVt5fu4vFPNxLfNoppd53C0ITYQIfVZAXqKqoj/jCEs/LKKp6at5lXv97BqJ5xvHTjaBu8LDS9AETinHYCWKqqd4pIN+C/qnoxEA9Md5c3B95V1XmBCrixlZRX8vinG3l32W5O7tuB564fZYOXmbBUVFbBw9OTmL46nbMHdeKf144krpWNth1IgeqDU+cPQ4BiaVRZB0u4993VfLczh8mn9OZ3F59Ay+Z2TjYUqWr/eh7fA1zs3k8BRjRmXMEi7UARd7+zinVpedx5Zj9+dcFAmlv/g7AlIlcAzwOdgNkiskZVLwxwWI0iZV8Bd729ii1Z+fzy/IHcc3Z/60QcBAJ1FVWdPwzhblnKfu55dzWFpRU8e91IJo1MCHRIxjSIL5OzeOCDNVRWKv+5+SQuPLGubhsmnKjqdGB6oONobPOS9vKrD9fRopnw5u3jOH1Ap0CHZFzBcBVV2FNVXv06hb/NS6ZX+1a8+3/jGWhT3ZswVFWlPPfFVp79fCuD4mN4+aaT6N2xdaDDMsbvKiqrePqzZP7zVQojejhdDRLirKtBMLEEp4Hll5Tz4IfrmLchg4uHdeFvPx5uV46YsHSgsIwHPljD4i37uHJ0An++fJhdOWLCUla+29VgRw43T+jFI5ecQGRzq+vBxhKcBpSckc+db69kd04Rj/zoBH5yWh8bBMuEpbWpudz9zir25ZfylyuGcf24HlbXTVj6bkcO97y7ivyScv557QiuGNU90CGZeliC00BmrE7noWnraRPVnPf+bwLj+thlsSb8qCrvfZfKY59soFNMJB/ddbJNjmnCkqry2pIdPDl3Mz3bt+Ktn4yzyTGDnCU4flZaUcmfZ2/izf/tYlyf9rxwwyg6x0QFOixj/K64rJJHZiTx8ao0zhjYiWevHUm71nZZrAk/+SXl/ObjdcxZn8GFJ8bz9NUjaGtdDYKeJTh+tCe3mLvfWcWa1FymnNGXX184yC6LNWFpZ3Yhd769kuTMfB44bwD3nTPA5pEyYWlLptPVYNf+In538WD+7/S+dvo1RFiC4ydLtmbzs/dXU1ZRxcs3jWbi0K6BDsmYBjF/Qwa//HAtzSKExMljOWtQ50CHZEyDmLkmnd9+vJ7Wkc15547xTOjbIdAhmaNgCY4fzFidzs+nrmFg5xj+fdNo+nZqE+iQjGkQby/dxSMzkhjePZaXbhxN93atAh2SMQ3i5cXb+evczYzt3Y4XbhhNfFvrahBqLME5Tt9sy+bBj9YyoU8HXps8hlYt7SM14Wn+hgx+PzOJcwd35sUbR9sM4CZsTVuVxl/nbubSEd145poRNgN4iLL/2nHYuOcgP31rJf06teE/t5xkyY0JWyt3HeC+91YzvHscL9xgyY0JX0u2ZvPrj9Zxct8O/P3q4ZbchDD7zx2j9NxibvN8R0xUcxJvG2s96k3YStlXwB1vLKdrbBSv3TrGBu8zYWvjnoPc+bZz0PryzSfZ4H0hzpocjkFeUTmTX/+OorJKPrrzFJsJ3IStffml3Jr4HREivHH7ODq0sZnATXiqPmhtE9kcz+1jiY22g9ZQZwnOUSqtqGTKWyvYtb8Iz+1jGdTF5pQy4amwtIKfvLGc7Pwy3psygV4dbE4pE55qDlpLK/nwrpPtoDVMWIJzFKqqlF9OXcuyHTk8e91ITunXMdAhGdMgKiqruPfdVSSl5/HqLWMY2SMu0CEZ0yBKKyr5v7dWsHN/IW/cbqMThxNLcI7Ck3M38em6vTx00WAmjUwIdDjGNAhV5ZEZSSxK3sdfrhjGuSfEBzokYxpEVZXyi6lr+c4OWsOSdTL20etLdvDq1zuYfEpvppzRN9DhmCAgIo+JSLqIrHFvF9ez3kQRSRaRbSLy28aO82g9/8U23l+eyn3n9OeG8T0DHY4xDebJuZuYbQetYctacHwwd/1e/jR7IxeeGM+jlwyxYbqNt3+q6t/rWygizYAXgfOBNGC5iHyiqhsbK8CjMXVFKs8s2MKPR3fnF+cPDHQ4xjSY6oPWW0/uZQetYcpacI5g+c4c7v9gDaN7tuPZ60bZfDvmaI0DtqlqiqqWAe8DkwIcU50Wb9nHQ9PWc/qAjvz1x8MskTdha47XQevvLz3R6nqYsgTnMLZl5XPHGyvoHhfNf28ZY4ObmbrcKyLrROR1EWlXx/IEINWrnOY+dggRmSIiK0Rkxb59+xoi1nolpedx19srGRQfw0s3jrbBzUzY+m5HDg98sIZRPeLsoDXMBXQvJiK/EhEVkaDr2ZV1sIRbX19Oi2YRvHH7ONq1bhnokEwAiMhCEUmq4zYJ+DfQDxgJ7AX+Udcm6nhM63otVX1FVceo6phOnTr56y0cUWpOEZMTl9OuVUsSbxtLjA1aaY6BiDwtIpvdhH+6iMQFOqbatmXl839vOgetr9061g5aw1zA+uCISA+cfgm7AxVDfQpKK7jNs5wDRWV8MOVkerS3CQWbKlU9z5f1RORV4NM6FqUBPbzK3YE9fgjNLw4UlnFr4neUVVTy/pTxNqGgOR4LgIdUtUJE/gY8BPwmwDHV+P6gVeygtYkIZAvOP4FfU8/RbKCUV1Zx19sr2ZyRz4s3jmZY99hAh2SClIh09SpeASTVsdpyYICI9BGRlsB1wCeNEd+RlJRXcsebK0g7UMx/bx1L/842aKU5dqo6X1Ur3OJSnGQ+KBSUVjA50TloTZw8zg5am4iAJDgichmQrqprfVi30folqCq//Xg9X2/N5skrhnH2oM4N+nom5D0lIutFZB1wNvBzABHpJiJzANwd/r3AZ8AmYKqqbghUwNUqq5QH3l/Dqt0H+Ne1IxnXp32gQzLh5XZgbn0LG3O/Xn3QmpxpB61NTYOdohKRhUCXOhY9DPwOuMCX7ajqK8ArAGPGjGnQ1p5/LtjCx6vSeOC8AVwztseRn2CaNFW9uZ7H9wAXe5XnAHMaK64jUVX+9OlG5m3I4NFLhnDxsK5HfpIxHH6/rqoz3XUeBiqAd+rbTmPt170PWp/68XA7aG1iGizBqa/vgogMA/oAa91L87oDq0RknKpmNFQ8R/Lust0898U2rh3Tg/vPHRCoMIxpcK9+nYLn253ccVoffnJan0CHY0LIkfqkicitwCXAuaoa8O4Hz9hBa5PW6J2MVXU9UJNGi8hOYIyqZjd2LNW+2JzJIzPWc9agTjxxxVAbE8GErU/W7uEvczbzo+Fd+d3FJwQ6HBNGRGQiTqfiM1W1KNDxvPfdbp63g9YmrckPdnGgsIxfTl3LCV3b8uINNv6HCV8ZeSU89PE6xvZuxz+uHkGEjf9h/OsFIAZY4E5d8nKgAtmRXcgfZm7gjIF20NqUBXyqBlXtHcjXf+qzzRwsqeC9a0bQOjLgH4cxDeZPszdSUaX84+qRNv6H8TtV7R/oGMDpd/OHTzbQsnkEf79quB20NmFN+j+/evcB3l+eym2n9GZwl7aBDseYBrNkazaz1+3l7rP607ODXSJrwte8pAy+2rKPn58/kM42rlOT1mQTnMoq5dGZSXSOieQBm1TQhLHSikp+PzOJ3h1a8dMzbVJBE74KSyt4/NONDO4Sw60n9wp0OCbAmmyC886yXSSlH+TRS4bQxk5NmTD23693kJJdyB8nDbVTUyasPffFVvbmlfDnK4bS3E5NNXlNsgbsyy/l6c+SOa1/R35kY4CYMJaaU8TzX2zloqFdOHNg481vZUxj25qZz2tf7+Dqk7pzUi8buNI00QTnyTmbKCmv5I+TTrTe9Sas/XHWRiJEePSSIYEOxZgGo6o8MiOJ1pHN+e1FgwMdjgkSTS7BWZayn2mr05lyRl/6dWoT6HCMaTCfb8pk4aZMfnbuALrFRQc6HGMazMw1e1i2I4cHLxxEhzaRgQ7HBIkmleCUV1bx6MwkEuKiufdsG/jJhK+S8koem7WB/p3bcPupNlqxCV8HS8p5YvYmhneP5fpxPQMdjgkiTap3reebnWzJLOCVm08iuqV1tjTh66VF20jNKea9/5tAy+ZN6jjGNDHPzN/C/sJSXp88hmY2eKXx0mT2fBl5Jfxr4RbOHdyZ84fEBzocYxrMjuxCXl6cwqSR3Ti5X4dAh2NMg9mwJ483/7eTG8f3ZHj3uECHY4JMk0lwqkdxfewy61hswlf1KK6RzSN42OaaMmGsqkp5dEYS7Vq15MELrGOxOVSTSHC+3rqP2ev2cs/Z/enR3kZxNeHLRnE1TcVHK9NYtTuX3140mNhWLQIdjglCYZ/gOKO4bqB3h1ZMOcNGcTXhq7C0gj/O2sgJXdtyi43iasLYgcIynpy7iTG92vHj0d0DHY4JUmHfyfjVr1LYkV3IG7ePs1FcjV+JyAfAILcYB+Sq6sg61tsJ5AOVQIWqjmmIeJ77fCsZB0t48cbRNoqrCWtPfZbMwZIK/nT5UCKsY7GpR1gnOKk5RbywaJuN4moahKpeW31fRP4B5B1m9bNVNbuhYtmSmc9rS3ZwzZjunNSrXUO9jDEBtyY1l/eX7+a2U/pwQlebJNnUL6wTHBvF1TQGcXqtXwOcE4jXV3U6W7aObM5vJlpnSxO+KquUR2asp1ObSH5+vo1lZg4vbNuxq0dxvd9GcTUN73QgU1W31rNcgfkislJEpvj7xatHcf31RBvF1YS3d91Jkh+5ZAgxUdax2BxeWLbgVI/iOqBzG24/zUZxNcdORBYCXepY9LCqznTvXw+8d5jNnKqqe0SkM7BARDar6ld1vNYUYApAz56+jchaPYrriO6xXDfWRnE14Su7wJkk+ZR+Hbh0uE2SbI4sLBMc71FcW1hnS3McVPW8wy0XkebAlcBJh9nGHvdvlohMB8YBhyQ4qvoK8ArAmDFj1Jf4bBRX01Q8OWczxeWVPD5pqI1lZnwSsF9/EblPRJJFZIOIPOWv7VaP4nq5jeJqGsd5wGZVTatroYi0FpGY6vvABUCSP17YRnE1wURE/iQi60RkjYjMF5Fu/tr2dzty+HhVGnec3pf+nW2SZOObgCQ4InI2MAkYrqonAn/3x3ZVld/PTCKyeQS/+5GN4moaxXXUOj0lIt1EZI5bjAeWiMha4DtgtqrOO94XrapSHrFRXE1weVpVh7tDJXwK/N4fGy2vrOLRGc4kyfed098fmzRNRKBOUd0F/FVVS8FpuvfHRucmZfD11mz+cOkQOsfYKK6m4anq5Doe2wNc7N5PAUb4+3U/XJnK6t25PH3VcBvF1QQFVT3oVWyN07n+uL3x7U6SM/N5+aaTaNUyLHtVmAYSqFNUA4HTRWSZiCwWkbH1rSgiU0RkhYis2LdvX70bLCyt4PFZGxnStS03T7BRXE34OlBYxl/nbmZsbxvF1QQXEfmziKQCN3KYFhxf9+sZeSX8c8EWzhrUiQtPtEmSzdFpsARHRBaKSFIdt0k4LUftgAnAg8BUqafXmKq+oqpjVHVMp071D9ZXPYrrny4faqO4mrBmo7iaQDnCfh1VfVhVewDvAPfWtx1f9+tPzN5IeZXyR5sk2RyDBmvvO9zVJyJyFzBNVRX4TkSqgI5A/an8YVSP4nrtmB42iqsJa9WjuN5+ah8Gd7FRXE3jOtJVhV7eBWYDfzjW1/pmWzafrtvLA+cNoFeH1se6GdOEBaqpYwbuqK8iMhBoCRzTMPbVo7i2iWrOby6yzpYmfHmP4vrAeTaKqwkuIuJdKS8DNh/rtkorKnl0ZhK9OrTizjP7HX9wpkkKVI+t14HXRSQJKANudVtzjtriLftYtiOHv1wxjPatW/o1SGOCyZz1e0lKP8hz14+yUVxNMPqriAwCqoBdwJ3HuqGPV6aTsq+QxNvG2iTJ5pgFJMFR1TLgJn9s68yBnXjt1jGcPaizPzZnTND60bCutGrZjHMGW103wUdVf+yvbV07tgddYiNtv26OS8hfcycinHuC9a434S8iwuq6aRqaRQjnDLa6bo6PXW5kjDHGmLBjCY4xxhhjwo4cY9/egBCRfTid1+rSkWO8EqsRWYz+cTwx9lLV+gfeCBKHqeuh8P+B0Igz3GO0ut44QiHOcI+xzroeUgnO4YjIClUdE+g4Dsdi9I9QiLGhhMp7D4U4LcbgFirvPRTibKox2ikqY4wxxoQdS3CMMcYYE3bCKcF5JdAB+MBi9I9QiLGhhMp7D4U4LcbgFirvPRTibJIxhk0fHGOMMcaYauHUgmOMMcYYA1iCY4wxxpgwFFIJjohMFJFkEdkmIr+tY7mIyHPu8nUiMrqR4+shIotEZJOIbBCR++tY5ywRyRORNe7t940Zo1ccO0VkvRvDijqWB/qzHOT1Ga0RkYMi8kCtdYLis2wIVtf9GqvV9SBmdd1vcQZ1PXdjaNy6rqohcQOaAduBvkBLYC0wpNY6FwNzAQEmAMsaOcauwGj3fgywpY4YzwI+DYLPcyfQ8TDLA/pZ1vG/z8AZzCnoPssGer9W1/0Xq9X1IL1ZXfdrnCFTz73+9w1a10OpBWccsE1VU9SZjfx9YFKtdSYBb6pjKRAnIl0bK0BV3auqq9z7+cAmIKGxXt/PAvpZ1nIusF1V6xvFOtxYXW9cVtcDx+p64wmmeg6NUNdDKcFJAFK9ymkcWsl8WadRiEhvYBSwrI7FJ4vIWhGZKyInNm5kNRSYLyIrRWRKHcuD5rMErgPeq2dZMHyW/mZ13b+srgcvq+v+E0r1HBqhrjc/1icGgNTxWO1r3H1Zp8GJSBvgY+ABVT1Ya/EqnCa5AhG5GJgBDGjkEAFOVdU9ItIZWCAim1X1K6/lwfJZtgQuAx6qY3GwfJb+ZnXdv6yuBy+r6/4TEvUcGq+uh1ILThrQw6vcHdhzDOs0KBFpgfMleEdVp9VerqoHVbXAvT8HaCEiHRszRve197h/s4DpOE3F3gL+WbouAlapambtBcHyWTYAq+t+ZHU9qFld95MQqufQSHU9lBKc5cAAEenjZn/XAZ/UWucT4Ba3t/gEIE9V9zZWgCIiwGvAJlV9pp51urjrISLjcP4H+xsrRvd1W4tITPV94AIgqdZqAf0svVxPPc2YwfBZNhCr635idT3oWV33T4yhVM+hkep6yJyiUtUKEbkX+Ayn9/XrqrpBRO50l78MzMHpKb4NKAJua+QwTwVuBtaLyBr3sd8BPb1ivAq4S0QqgGLgOlVt7GbCeGC6W4eaA++q6rwg+ywRkVbA+cBPvR7zjjEYPku/s7ruV1bXg5jVdb8JiXoOjVvXbaoGY4wxxoSdUDpFZYwxxhjjE0twjDHGGBN2LMExxhhjTNixBMcYY4wxYccSHGOMMcaEHUtwjDHGGBN2LMExxhhjTNj5f4KZvGuIEuqrAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tfp = (irf['C'] + ss['C'])/(irf['L'] + ss['L'])\n",
    "fig, axes = plt.subplots(2, 3, figsize=(8, 4))\n",
    "ax = axes.flatten()\n",
    "\n",
    "ax[0].plot(100*((irf['N'][1:10] + ss['N'])/ss['N']-1))\n",
    "ax[0].axhline(0, color='gray', linestyle=':')\n",
    "ax[0].set_title('Mass of establishments')\n",
    "ax[0].set_ylabel('pct deviation from ss')\n",
    "\n",
    "ax[1].plot(100*((irf['mu'][1:10] + ss['mu'])/ss['mu']-1))\n",
    "ax[1].axhline(0, color='gray', linestyle=':')\n",
    "ax[1].set_title('Markup')\n",
    "\n",
    "ax[2].plot(100*((tfp[1:10])/tfp[-1]-1))\n",
    "ax[2].axhline(0, color='gray', linestyle=':')\n",
    "ax[2].set_title('TFP')\n",
    "\n",
    "ax[3].plot(100*((irf['L'][1:10] + ss['L'])/ss['L']-1))\n",
    "ax[3].axhline(0, color='gray', linestyle=':')\n",
    "ax[3].set_title('Employment')\n",
    "\n",
    "ax[4].plot(100*((irf['C'][1:10] + ss['C'])/ss['C']-1))\n",
    "ax[4].axhline(0, color='gray', linestyle=':')\n",
    "ax[4].set_title('Consumption')\n",
    "\n",
    "ax[5].plot(100*((irf['w'][1:10] + ss['w'])/ss['w']-1))\n",
    "ax[5].axhline(0, color='gray', linestyle=':')\n",
    "ax[5].set_title('Wage')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "232e60a1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3e395ab550>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAUfElEQVR4nO3db4xcV3nH8e/j9a6dODRxcP5sbCc2wWpr0UDDylAVoVZJqGNaHFohOVKL1SJFVE3/vEDFKBUFVZXSohapagR1WyS3oo1SWogLRsFYrfKiBWJDEuIY400gxLGxnTRxYjvenT9PX8zd3dn1rHftsT27nO9HWs3cO+fuPDc3Pr8999w7E5mJJKlcC3pdgCSptwwCSSqcQSBJhTMIJKlwBoEkFW5hrws4H8uWLctVq1b1ugxJmlf27NnzYmZeM3X9vAyCVatWsXv37l6XIUnzSkQ812m9p4YkqXAGgSQVziCQpMIZBJJUOINAkgp3QYIgItZHxP6IGI6ILR1ej4j4m+r1JyPi1tluK0m6uLoOgojoAx4A7gTWAndHxNopze4E1lQ/9wCfOYdtJUkX0YW4j2AdMJyZzwJExIPARuDptjYbgX/K1mdefyMiroqIQWDVLLadU06M1Dn48ikOHz/NqZEGp0brvF5rcLrWoNGEJMmEzOoRqseJ9ZJ0vt5/6wpWL1tyQX/nhQiC5cDzbcsHgXfMos3yWW4LQETcQ2s0wY033thdxefo6UOv8m97nud/hl9i/5HXuv59ERegKElFuvWmpXMyCDp1a1P/7J2uzWy2ba3M3ApsBRgaGrokf1a/drrGJ//zab6w5yCLFi5g3eqr+bW3DrJq2RIGr1zMFYv6uXygj8sG+ljc30dfBBGtjj6ongMRUT22nkvSXHIhguAgsLJteQVwaJZtBmaxbU8cP1Xjtz73TfYeepXf/aWb+fC7b+bKy/t7XZYkXXAX4qqhx4A1EbE6IgaATcD2KW22Ax+srh56J3A8Mw/PcttLLjP5yBeeYN/hV/m733w7H13/M4aApJ9YXY8IMrMeEfcCjwB9wOcyc29EfLh6/bPADmADMAycAn77bNt2W1O3vvT4C+x8+gh/8t6f5fa11/W6HEm6qGI+XsUyNDSUF+vTR+uNJnd8+lEu6+/jy7//LhYs8Jy+pJ8MEbEnM4emrvfO4im++tSP+cGLJ/mD29YYApKKYBBM8aXvvMANVy7mPZ4SklQIg6DN8ddrPHrgGO+9ZdDRgKRiGARtdu07Qq2RvPeWG3pdiiRdMgZBm/995iWWXt7PW1dc2etSJOmSMQja7HnuZd5+01Lv/pVUFIOg8tKJEZ598SRvv+nqXpciSZeUQVDZ89zLAAytWtrjSiTp0jIIKk+9cJwFAT+33PkBSWUxCCrDx05w0xuXsLi/r9elSNIlZRBUho+e4OZrLuxnfEvSfGAQ0Pp8oR++eIqbr72i16VI0iVnEADPv/w6o40mN19jEEgqj0EAPHP0BABvdkQgqUAGAfDMsVYQ3LzMIJBUHoMAOPTK67xh8UK/hUxSkQwC4MevnmbwysW9LkOSesIgAH58/DTX/ZRBIKlMBgFw+LgjAknlKj4Iao0mx06McP2Vl/W6FEnqieKD4NhrI2TiiEBSsYoPgsPHTwNwvXMEkgpVfBAcebUKAkcEkgpVfBCMjQg8NSSpVMUHwbHXRujvC668zJvJJJWp+CB45dQoV10+4PcUSypW8UHw8qlRlvrREpIKZhCcqnHV5QO9LkOSeqb4IHjFEYGkwhUfBC+fqrHUEYGkghUdBJk5PlksSaUqOghOjjaoNdJTQ5KKVnQQvHxyFMBTQ5KKVnQQvHKqBsBVjggkFazoIHj5VDUiWOKIQFK5DAJwjkBS0boKgoi4OiJ2RsSB6nHpNO3WR8T+iBiOiC1t6z8QEXsjohkRQ93Ucj4mTg05IpBUrm5HBFuAXZm5BthVLU8SEX3AA8CdwFrg7ohYW738FPDrwKNd1nFexoPAD5yTVLBug2AjsK16vg24q0ObdcBwZj6bmaPAg9V2ZOa+zNzfZQ3n7bXTNS4f6GNhX9FnyCQVrtse8LrMPAxQPV7boc1y4Pm25YPVunMSEfdExO6I2H3s2LHzKnaqk6N1lixaeEF+lyTNVzP2ghHxdeD6Di/dN8v36PT5zjnLbSc2yNwKbAUYGho65+07OTHS4AqDQFLhZuwFM/P26V6LiCMRMZiZhyNiEDjaodlBYGXb8grg0DlXehGcHKmzZFFfr8uQpJ7q9tTQdmBz9Xwz8HCHNo8BayJidUQMAJuq7XruxEidJQOOCCSVrdsguB+4IyIOAHdUy0TEDRGxAyAz68C9wCPAPuChzNxbtXt/RBwEfgH4SkQ80mU95+TkSN1TQ5KK11UvmJkvAbd1WH8I2NC2vAPY0aHdF4EvdlNDN1qnhgwCSWUr+rrJEyMNg0BS8YoOgtapISeLJZWt2CBoNJPXa44IJKnYIDg5WgdwslhS8YoNghOnW0HgiEBS6YoNgpMjBoEkQcFBcGJk7NSQk8WSylZsEJwcaQB4Z7Gk4hUbBCc8NSRJQMFBcHLEq4YkCUoOglFHBJIEBQfBCUcEkgQUHAQnR+osCFjcX+x/AkkCCg6C10ebXNbfR0SnL1CTpHIUGwSn6w0W93sPgSSVGwQ1g0CSoOAgGKk3WeT8gCQVHAS1BosXOiKQpGKD4HSt6RVDkkTRQeAcgSRByUHgVUOSBJQcBJ4akiSg6CBwsliSoOggaLLIU0OSVG4QjNQaLFpY7O5L0rhie0IniyWppcggaDSTWiOdLJYkCg2C07XW9xU7IpCk0oPAOQJJKjQI6k3AEYEkQaFBMOKpIUkaV2QQnK6NjQiK3H1JmqTInvB0vTUi8IYySSo1CMYniw0CSSoyCEY8NSRJ44rsCb2PQJImdBUEEXF1ROyMiAPV49Jp2q2PiP0RMRwRW9rWfyoivhcRT0bEFyPiqm7qma2xOQKDQJK6HxFsAXZl5hpgV7U8SUT0AQ8AdwJrgbsjYm318k7gLZl5C/B94GNd1jMrXjUkSRO67Qk3Atuq59uAuzq0WQcMZ+azmTkKPFhtR2Z+LTPrVbtvACu6rGdWnCyWpAndBsF1mXkYoHq8tkOb5cDzbcsHq3VT/Q7w1S7rmZWJEYFBIEkLZ2oQEV8Hru/w0n2zfI/osC6nvMd9QB34/FnquAe4B+DGG2+c5Vt3NjYi8PsIJGkWQZCZt0/3WkQciYjBzDwcEYPA0Q7NDgIr25ZXAIfafsdm4FeB2zIzmUZmbgW2AgwNDU3bbjZO1xsM9C1gwYJOGSVJZen2T+LtwObq+Wbg4Q5tHgPWRMTqiBgANlXbERHrgY8C78vMU13WMmuj9SYDjgYkCeg+CO4H7oiIA8Ad1TIRcUNE7ACoJoPvBR4B9gEPZebeavu/Bd4A7IyIxyPis13WMyu1hkEgSWNmPDV0Npn5EnBbh/WHgA1tyzuAHR3avbmb9z9ftXoy0GcQSBIUemdxrdGkf6HzA5IEhQbBaKNJvyMCSQIKDYJao+mpIUmqFNkb1hrpiECSKkX2hrVGk/4+5wgkCQoNgtG6cwSSNKbI3nDU+wgkaVyRvWHNq4YkaVyRvWGtns4RSFKlzCBoNBnwuwgkCSg0CEa9akiSxhUZBN5QJkkTiuwNvaFMkiYU2RvWvI9AksYV2RuO+umjkjSuyCBwjkCSJhTXG9YbTZqJp4YkqVJcb1hrtL733iCQpJbiesPRRhPA+wgkqVJcENSqIFjkh85JElBwEHhqSJJaiusNa3XnCCSpXXG94fgcgaeGJAkoMAjGTg0NOFksSUDBQeCpIUlqKa43NAgkabLiesNRJ4slaZLiesOxyeIBP3ROkoACg6BW99SQJLUrrjd0jkCSJiuuN5w4NVTcrktSR8X1hmOfPur3EUhSS3G9oaeGJGmy4nrDmh9DLUmTFBcEo3U/a0iS2hXXGzpHIEmTddUbRsTVEbEzIg5Uj0unabc+IvZHxHBEbGlb/2cR8WREPB4RX4uIG7qpZzacI5CkybrtDbcAuzJzDbCrWp4kIvqAB4A7gbXA3RGxtnr5U5l5S2a+Dfgy8PEu65nRaL3JgoC+Bc4RSBJ0HwQbgW3V823AXR3arAOGM/PZzBwFHqy2IzNfbWu3BMgu65lRrdF0NCBJbRZ2uf11mXkYIDMPR8S1HdosB55vWz4IvGNsISL+HPggcBz45S7rmdFoo+n8gCS1mbFHjIivR8RTHX42zvI9Op2DGf/LPzPvy8yVwOeBe89Sxz0RsTsidh87dmyWb32mRjNZ6KWjkjRuxhFBZt4+3WsRcSQiBqvRwCBwtEOzg8DKtuUVwKEO7f4F+Arwp9PUsRXYCjA0NHTep5DqzaRvgSMCSRrTbY+4HdhcPd8MPNyhzWPAmohYHREDwKZqOyJiTVu79wHf67KeGdUbTRY6USxJ47qdI7gfeCgiPgT8CPgAQHUZ6D9k5obMrEfEvcAjQB/wuczcO7Z9RPw00ASeAz7cZT0zao0IDAJJGtNVEGTmS8BtHdYfAja0Le8AdnRo9xvdvP/5aDTTj5eQpDbFnSx3RCBJkxUXBI1GstDJYkkaV1yPWG82HRFIUpsCg8D7CCSpXXFB0Giml49KUpvigqDuHIEkTVJcj9jwqiFJmqS4IKg3m84RSFKbAoPAEYEktSsvCJwjkKRJiusRvWpIkiYrLgjqzSZ9zhFI0rjigsARgSRNVlwQ1BpOFktSu+KCoNFM+p0slqRxxfWI9WY6RyBJbYoLgkbTr6qUpHbFBYE3lEnSZOUFQcOrhiSpXXFB0PrQueJ2W5KmVVyPWG82/fJ6SWpTVBA0m0kzcY5AktoUFQSNTADnCCSpTVFBUG+0gsA5AkmaUFSPWG82AUcEktSuqCBoNKtTQ04WS9K4ooKg3nSOQJKmKioIxkYEzhFI0oSiekRHBJJ0prKCoNGaLPY+AkmaUFYQOFksSWcoKgjGrxpyjkCSxhXVI07cUOaIQJLGFBUEDSeLJekMRQVBrbqz2K+qlKQJRQWBIwJJOlNRQTA2R+BksSRN6KpHjIirI2JnRByoHpdO0259ROyPiOGI2NLh9Y9EREbEsm7qmYmfNSRJZ+r2T+MtwK7MXAPsqpYniYg+4AHgTmAtcHdErG17fSVwB/CjLmuZ0dinj3rVkCRN6DYINgLbqufbgLs6tFkHDGfms5k5CjxYbTfm08AfA9llLTOaODVkEEjSmG6D4LrMPAxQPV7boc1y4Pm25YPVOiLifcALmfnETG8UEfdExO6I2H3s2LHzKrbe9D4CSZpq4UwNIuLrwPUdXrpvlu/RqdfNiLi8+h3vmc0vycytwFaAoaGh8xo9jM0R9Pc5WSxJY2YMgsy8fbrXIuJIRAxm5uGIGASOdmh2EFjZtrwCOATcDKwGnoiIsfXfjoh1mfnjc9iHWXOOQJLO1O2fxtuBzdXzzcDDHdo8BqyJiNURMQBsArZn5ncz89rMXJWZq2gFxq0XKwTA+wgkqZNug+B+4I6IOEDryp/7ASLihojYAZCZdeBe4BFgH/BQZu7t8n3Pi3MEknSmGU8NnU1mvgTc1mH9IWBD2/IOYMcMv2tVN7XMhjeUSdKZiuoRG9UcgTeUSdKEooLAr6qUpDMVFQQN5wgk6QxFBUHdbyiTpDMU1SP65fWSdKaygsA5Akk6Q1FB0GgmCwIWGASSNK6oIKg30/kBSZqiqF6x0UznByRpiqKCoNZoOj8gSVMUFQSNZtLnXcWSNElXnzU036wd/ClO1xq9LkOS5pSigmDTuhvZtO7GXpchSXNKUaeGJElnMggkqXAGgSQVziCQpMIZBJJUOINAkgpnEEhS4QwCSSpcZGavazhnEXEMeO48N18GvHgBy+kl92Vucl/mJvcFbsrMa6aunJdB0I2I2J2ZQ72u40JwX+Ym92Vucl+m56khSSqcQSBJhSsxCLb2uoALyH2Zm9yXucl9mUZxcwSSpMlKHBFIktoYBJJUuKKCICLWR8T+iBiOiC29rudcRcQPI+K7EfF4ROyu1l0dETsj4kD1uLTXdXYSEZ+LiKMR8VTbumlrj4iPVcdpf0T8Sm+qPtM0+/GJiHihOi6PR8SGttfm5H4ARMTKiPiviNgXEXsj4g+r9fPxuEy3L/Pu2ETE4oj4VkQ8Ue3LJ6v1F++4ZGYRP0Af8AzwJmAAeAJY2+u6znEffggsm7LuL4Et1fMtwF/0us5pan83cCvw1Ey1A2ur47MIWF0dt75e78NZ9uMTwEc6tJ2z+1HVNwjcWj1/A/D9qub5eFym25d5d2yAAK6onvcD3wTeeTGPS0kjgnXAcGY+m5mjwIPAxh7XdCFsBLZVz7cBd/WulOll5qPA/01ZPV3tG4EHM3MkM38ADNM6fj03zX5MZ87uB0BmHs7Mb1fPXwP2AcuZn8dlun2Zzlzel8zME9Vif/WTXMTjUlIQLAeeb1s+yNn/R5mLEvhaROyJiHuqdddl5mFo/WMAru1Zdeduutrn47G6NyKerE4djQ3Z581+RMQq4Odp/fU5r4/LlH2BeXhsIqIvIh4HjgI7M/OiHpeSgiA6rJtv187+YmbeCtwJ/F5EvLvXBV0k8+1YfQa4GXgbcBj4q2r9vNiPiLgC+HfgjzLz1bM17bBuTu1Ph32Zl8cmMxuZ+TZgBbAuIt5yluZd70tJQXAQWNm2vAI41KNazktmHqoejwJfpDX8OxIRgwDV49HeVXjOpqt9Xh2rzDxS/cNtAn/PxLB8zu9HRPTT6jg/n5n/Ua2el8el077M52MDkJmvAP8NrOciHpeSguAxYE1ErI6IAWATsL3HNc1aRCyJiDeMPQfeAzxFax82V802Aw/3psLzMl3t24FNEbEoIlYDa4Bv9aC+WRn7x1l5P63jAnN8PyIigH8E9mXmX7e9NO+Oy3T7Mh+PTURcExFXVc8vA24HvsfFPC69niG/xLPxG2hdTfAMcF+v6znH2t9E68qAJ4C9Y/UDbwR2AQeqx6t7Xes09f8rraF5jdZfMB86W+3AfdVx2g/c2ev6Z9iPfwa+CzxZ/aMcnOv7UdX2LlqnEJ4EHq9+NszT4zLdvsy7YwPcAnynqvkp4OPV+ot2XPyICUkqXEmnhiRJHRgEklQ4g0CSCmcQSFLhDAJJKpxBIEmFMwgkqXD/D0ZdKuhcYiLMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(irf['L'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec0593fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3e3952a160>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD6CAYAAACxrrxPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAf1ElEQVR4nO3deXhVhZ3/8fc3NysQ1gTQsAQUUUTcAqIVl6m2gI6OtVXR2mpLqU7tdKYzv6p1WjtjbV3q1GlrddCiXazaqbZaQa12GaiKEBVZZA8Bwhr2QAJZ7vf3R24wxpAc4N6cu3xez3Of3HPPyb0f8sCHk3PPuV9zd0REJPVlhR1ARETiQ4UuIpImVOgiImlChS4ikiZU6CIiaUKFLiKSJjotdDObYWZbzWxxJ9uNNbMmM/t0/OKJiEhQ1tl56GZ2HrAX+IW7jz7ENhHgVWA/MMPdf9vZCxcVFXlpaelhBxYRyWRvv/32Nncvbm9ddmff7O6zzay0k82+CjwLjA0aqrS0lPLy8qCbi4gIYGZrD7XuqI+hm1kJcAXwyNE+l4iIHLl4vCn6IHCruzd1tqGZTTOzcjMrr66ujsNLi4hIi04PuQRQBjxtZgBFwGQza3T337fd0N2nA9MBysrK9CEyIiJxdNSF7u7DWu6b2RPAi+2VuYiIJFanhW5mTwEXAEVmVgXcCeQAuLuOm4uIJIkgZ7lMCfpk7n7DUaUREZEjpitFRUTShApdRKQLPfjaCuat2ZGQ51ahi4h0kZVbanjwtZXMrdiekOdXoYuIdJHH5qwhPyeLz44fmpDnV6GLiHSB6poD/O7dDVx5xiD6ds9NyGuo0EVEusAv36ykIRrli+cO63zjI6RCFxFJsLr6Jn45dy0XnTSA4cU9EvY6KnQRkQT77TtV7Kxt4EsThif0dVToIiIJFI06M/62hlMH9WJsaZ+EvpYKXUQkgV5buoU12/YxdcJwYh9imDAqdBGRBHpszhpKehcwafTAhL+WCl1EJEEWrN/FvModfOHcYWRHEl+3KnQRkQR5dE4FhfnZXD12cJe8ngpdRCQB1u+o5aVFm7h23BB65MVjllDnVOgiIgnw+OuVZJlxw8dKu+w1VegiInG2u66BZ+av49Ixx3BMr4Iue10VuohInD09bx376puYmuALidpSoYuIxFF9Y5THX6/knOP6MbqkV5e+tgpdRCSOZi7ayOY9+xN+mX97VOgiInHi7jw6ew3H9+/B+ScUd/nrq9BFROLkzdXbeX/THqaeO4ysrMRe5t+eTgvdzGaY2VYzW3yI9deZ2cLY7Q0zOzX+MUVEkt+jcyoo6pHLP5xeEsrrB9lDfwKY2MH6NcD57j4GuAuYHodcIiIpZeWWGv6yvJrPnV1Kfk4klAydXr7k7rPNrLSD9W+0WpwLDIpDLhGRlJLoeaFBxPsY+heBlw610symmVm5mZVXV1fH+aVFRMLRFfNCg4hboZvZhTQX+q2H2sbdp7t7mbuXFRd3/TvAIiKJ0BXzQoOIyyfGmNkY4DFgkrtvj8dzioikgpZ5oR8/MbHzQoM46j10MxsCPAdc7+4rjj6SiEjqeDY2L3TaeV1/IVFbne6hm9lTwAVAkZlVAXcCOQDu/gjwbaAf8NPYeKVGdy9LVGARkWQRjTo/66J5oUEEOctlSifrpwJT45ZIRCRFtMwL/fGU0xM+LzQIXSkqInKEunJeaBAqdBGRI/BebF7ojR8r7ZJ5oUEkRwoRkRTTMi/0mnFDwo5ykApdROQwrd9Ry6wunhcahApdROQwhTEvNAgVuojIYQhrXmgQKnQRkcMQ1rzQIFToIiIBNTRFeeKNcOaFBqFCFxEJaObCTWzaHc680CBU6CIiAbg702dXhDYvNAgVuohIAGHPCw1ChS4iEkDY80KDUKGLiHRi1dbmeaHXjw9vXmgQKnQRkU48NmcNedlZXH92ePNCg1Chi4h0oLrmAM+9s4FPnxnuvNAgVOgiIh1IlnmhQajQRUQOIZnmhQahQhcROYSWeaFfmpD8e+egQhcRaVc06syIzQsdN6xv2HECUaGLiLTjT8u2UrFtH1MnDE+KeaFBdFroZjbDzLaa2eJDrDcz+5GZrTKzhWZ2Rvxjioh0rUdnVyTVvNAgguyhPwFM7GD9JGBE7DYNePjoY4mIhCcZ54UG0WlSd58N7Ohgk8uBX3izuUBvMzsmXgFFRLrao3MqKMzL5uqxg8OOclji8V9PCbC+1XJV7DERkZRTtbOWlxZv5tqzhlCYnxN2nMMSj0Jv790Cb3dDs2lmVm5m5dXV1XF4aRGR+Hr89UoMkm5eaBDxKPQqoPXvJYOAje1t6O7T3b3M3cuKi5Pz84RFJHPtrmvg6XnJOS80iHgU+gvA52Jnu4wHdrv7pjg8r4hIl0rmeaFBZHe2gZk9BVwAFJlZFXAnkAPg7o8As4DJwCqgFrgxUWFFRBKlZV7o2cOTc15oEJ0WurtP6WS9A1+JWyIRkRC0zAu9+4rRYUc5YqlzgqWISIK4O4/OaZ4XesEJ/cOOc8RU6CKS8d5cvZ0lG5N7XmgQKnQRyXipMC80CBW6iGS0VJkXGoQKXUQyWsu80M+OHxJ2lKOmQheRjFVdc4Dn3m2eF9qvR17YcY6aCl1EMtYv566loSk15oUGoUIXkYxUV9/EL9+sTJl5oUGo0EUkI6XavNAgVOgiknFa5oWOSaF5oUGo0EUk47y2dAsV2/bxpRSaFxqECl1EMkpT1PmvV1cwpG+3lJoXGoQKXUQyyv+Wr2fZ5hpunXhiSs0LDSK9/jQiIh3Yd6CRB15dwZlD+zD5lPTaOwcVuohkkP+ZXUF1zQHuuOSktDp23kKFLiIZYfPu/UyfvZpLxxzDGUP6hB0nIVToIpIRfvDH5USjcOvEE8OOkjAqdBFJe0s27ubZd6q48WOlDO7bLew4CaNCF5G05u7cPXMpvQty+McLjw87TkKp0EUkrf152VbeWL2dr318BL0KcsKOk1CBCt3MJprZcjNbZWa3tbO+l5n9wczeM7MlZnZj/KOKiByexqYo35u1lOFF3blu/NCw4yRcp4VuZhHgIWASMAqYYmaj2mz2FeB9dz8VuAB4wMxy45xVROSwPDV/Paur93HbpBPJSbOLiNoT5E84Dljl7hXuXg88DVzeZhsHCq35xM4ewA6gMa5JRUQOQ83+Bh58dQVnDevLxaMGhB2nS2QH2KYEWN9quQo4q802PwFeADYChcDV7h6NS0IRkSPw07+uZvu+ep64ZFRaXkTUniB76O39JLzN8ieBBcCxwGnAT8ys50eeyGyamZWbWXl1dfVhRhURCaZqZy0/+9sarji9hFMG9Qo7TpcJUuhVwOBWy4No3hNv7UbgOW+2ClgDfOTsfXef7u5l7l5WXFx8pJlFRDp0/yvLMeD/fXJk2FG6VJBCnw+MMLNhsTc6r6H58Epr64CPA5jZAGAkUBHPoCIiQSxYv4vnF2xk6oRhHNu7IOw4XarTY+ju3mhmtwCvABFghrsvMbObYusfAe4CnjCzRTQfornV3bclMLeIyEe4O9+buZSiHrncfEF6X0TUniBviuLus4BZbR57pNX9jcAn4htNROTwvLJkC/Mqd3D3FaPpkReo3tJK+p+YKSIZob4xyj0vLWVE/x5cXTa4829IQyp0EUkLv5q7lsrttXxz8klpN4koqMz8U4tIWtld28CP/rySc48v4oKRmXsGnQpdRFLej/+8kt11DXxzcnpOIgpKhS4iKW3d9lp+/mYlnzlzEKOO/cj1jBlFhS4iKe3el5eRnZXFv34isy4iao8KXURS1ttrdzBz0Sa+fP5wBvTMDztO6FToIpKS3J3vzlxK/8I8pp03POw4SUGFLiIp6cWFm3h33S7+7RMj6ZabeRcRtUeFLiIpZ39DE/e+vIwTBxZy5ZmDwo6TNFToIpJyfvFmJVU76/j3S0YRycrc0xTbUqGLSErZsa+eH/95FReOLObcEUVhx0kqKnQRSSk/+tNKauub+Obkk8KOknRU6CKSMlZX7+VXc9dyzdjBjBhQGHacpKNCF5GUEI06tz+7iILcCP980Qlhx0lKKnQRSQm/eLOSeZU7uPPvT6a4MC/sOElJhS4iSW/t9n3c+/JyLhxZzJVnlIQdJ2mp0EUkqUWjzjd+u5DsiPH9T43J6E9T7IwKXUSS2q/eWstba3bwrUtHMbCXPq+lIyp0EUla67bXcs9Lyzj/hGI+oytCO6VCF5GkFI0633j2PSJmfP9Tp+hQSwCBCt3MJprZcjNbZWa3HWKbC8xsgZktMbP/i29MEck0T85bx9yKHfz7pSdxbO+CsOOkhE4/oszMIsBDwMVAFTDfzF5w9/dbbdMb+Ckw0d3XmVn/BOUVkQywfkct35+1lAkjiriqbHDYcVJGkD30ccAqd69w93rgaeDyNttcCzzn7usA3H1rfGOKSKZoOasly4x7rtRZLYcjSKGXAOtbLVfFHmvtBKCPmf3VzN42s8+190RmNs3Mys2svLq6+sgSi0ha+/W8dbxZsZ07LjmJEh1qOSxBCr29/x69zXI2cCZwCfBJ4Ftm9pFrc919uruXuXtZcXHxYYcVkfTWcqjl3OOLuGasDrUcriBjPqqA1j/ZQcDGdrbZ5u77gH1mNhs4FVgRl5QikvbcndueWwjAPVfqrJYjEWQPfT4wwsyGmVkucA3wQpttngcmmFm2mXUDzgKWxjeqiKSzp+at5/VV2/nmJScxqE+3sOOkpE730N290cxuAV4BIsAMd19iZjfF1j/i7kvN7GVgIRAFHnP3xYkMLiLpo2pnLXfPfJ9zjuvHteOGhB0nZQWarOrus4BZbR57pM3y/cD98YsmIpnA3bn9uUU4cK/OajkqulJUREL1zPz1zFm5jdsnn8TgvjrUcjRU6CISmg276vjuzKWcPbwf1+lQy1FToYtIKFoOtUTdue/TY8jK0qGWo6VCF5FQTJ9dwewV1dw+6UQdaokTFbqIdLm3KrZz3yvLueSUY/js+KFhx0kbKnQR6VJb9+znlqfeZWjfbrqAKM4CnbYoIhIPjU1RbnnqXfbub+RXXzyLwvycsCOlFRW6iHSZ+19Zzrw1O/jh1acycmBh2HHSjg65iEiXeGXJZv5ndgWfHT+EK07XOLlEUKGLSMJVbtvHv/3mPU4d1ItvXToq7DhpS4UuIgm1v6GJm598h0jEeOi6M8jLjoQdKW3pGLqIJNS3fr+YZZv3MOOGsfoUxQTTHrqIJMwz89fxv29X8dULj+fCkRo1nGgqdBFJiMUbdvOt55cwYUQRX7voIwPMJAFU6CISd7trG7j5ybfp1z2X/77mdCL6nJYuoWPoIhJXTVHn679ZwObd+3nmy2fTt3tu2JEyhvbQRSRu3J3//MMS/rRsK9++dBRnDOkTdqSMokIXkbj52d/W8PM31/KlCcO4/uzSsONkHBW6iMTFzIWb+O7MpUw+ZSC3Tzop7DgZSYUuIketvHIH//KbBZQN7cN/XXWahlWEJFChm9lEM1tuZqvM7LYOthtrZk1m9un4RRSRZLa6ei9Tf1FOSe8CHv1cGfk5uhI0LJ0WuplFgIeAScAoYIqZfeTDGGLb3Qu8Eu+QIpKcqmsOcMPj84iY8cSNY+mjM1pCFWQPfRywyt0r3L0eeBq4vJ3tvgo8C2yNYz4RSVK19Y1M/fl8qmsO8LMbxjK0X/ewI2W8IIVeAqxvtVwVe+wgMysBrgAeiV80EUlWTVHnn55awKINu/nxlDM4bXDvsCMJwQq9vXc3vM3yg8Ct7t7U4ROZTTOzcjMrr66uDhhRRJKJu/OdF5bw2tItfOeyk7l41ICwI0lMkCtFq4DBrZYHARvbbFMGPB2bDVgETDazRnf/feuN3H06MB2grKys7X8KIpLk3J17XlrGL+euZdp5w/mczjVPKkEKfT4wwsyGARuAa4BrW2/g7sNa7pvZE8CLbctcRFKbu/ODPy7nf2ZXcP34odw+6cSwI0kbnRa6uzea2S00n70SAWa4+xIzuym2XsfNRTLAg6+t5KG/rGbKuCH8x2UnE/uNXJJIoA/ncvdZwKw2j7Vb5O5+w9HHEpFk8uM/reS//7SSq8oGcfc/jNaFQ0lKV4qKSIce/utqHnh1BZ86o4Tvf2qMyjyJqdBF5JAem1PBvS8v4/LTjuX+T5+qzzVPcip0EWnXz/62hu/OXMolY47hgc+ozFOBBlyIyIe4Oz98dQU/+vMqJo0eyINXn0Z2RPt+qUCFLiIHNUWdO19YzK/mruPqssHcfcVolXkKUaGLCAD1jVG+/psFvLhwEzedfxy3ThypUxNTjApdRKitb+TLv3ybOSu38c3JJzLtvOPCjiRHQIUukuF21dZz4xPzWVi1m/s+PYarygZ3/k2SlFToIhmsonovU39eTtWuOh6+7gw+cfLAsCPJUVChi2SoN1Zt4+Yn3yGSZfx66lmUlfYNO5IcJRW6SAZ68q213Pn8EoYXd+dnnx/L4L7dwo4kcaBCF8kgjU1R7p61lMdfr+TCkcX8aMrpFObnhB1L4kSFLpIhdu6r52vPLGD2imq+8LFh3HHJSbr6M82o0EUywHvrd/GPT75Ddc0Bvv+pU5gybkjYkSQBVOgiaczd+dVb67jrD+9TXJjHb28+mzGDeocdSxJEhS6SpmrrG7njd4v53bsbuGBkMQ9efRq9u+WGHUsSSIUukoYWVe3ma0+/y5rt+/jXi0/gKxcer88xzwAqdJE0Eo060+dU8MAfl9Ovex5PTj2Lc44rCjuWdBEVukia2LS7jq8/8x5vVmxn8ikD+d4Vp+gQS4ZRoYukOHfnf8uruOvF92ly574rx/CZskH6pMQMFOiDjs1sopktN7NVZnZbO+uvM7OFsdsbZnZq/KOKSFsbd9Xx+cfn841nF3LSsT2Z9U8TuGrsYJV5hup0D93MIsBDwMVAFTDfzF5w9/dbbbYGON/dd5rZJGA6cFYiAotI8yCKX7+1lntfXk7Unf+47GSuHz9Ub3xmuCCHXMYBq9y9AsDMngYuBw4Wuru/0Wr7ucCgeIYUkQ+8vXYn335+MUs27uFjx/fj+1eMYUg/fRaLBCv0EmB9q+UqOt77/iLw0tGEEpGP2r73APe+vIzflFcxoGceP7n2dC455RgdXpGDghR6e39bvN0NzS6kudDPPcT6acA0gCFDdOmxSBBNUefX89bxg1eWs+9AI18+bzhf/fgIeuTpnAb5sCB/I6qA1iNMBgEb225kZmOAx4BJ7r69vSdy9+k0H1+nrKys3f8UROQDcyu2c/fMpSzasJuzh/fjPy8/mREDCsOOJUkqSKHPB0aY2TBgA3ANcG3rDcxsCPAccL27r4h7SpEMs2Tjbu57eTn/t6KagT3z+dGU0/n7MTq8Ih3rtNDdvdHMbgFeASLADHdfYmY3xdY/Anwb6Af8NPYXrtHdyxIXWyQ9rdteywOvLuf5BRvpVZDD7ZNO5PPnlJKfEwk7mqQAcw/nyEdZWZmXl5eH8toiyWbT7joe/utqfv3WOrIjxo0fG8ZN5x9HrwINn5APM7O3D7XDrHdVREK0dvs+Hv7rap59p4qow9VjB/O1j49gQM/8sKNJClKhi4RgxZYaHvrLKv7w3kayI1lcPXYwXz7vOM32lKOiQhfpIu7OnJXbePz1NfxleTXdciNMnTCcqecOo7/2yCUOVOgiCVZX38Rz71bxxOuVrNy6l6IeefzzRSP4/Nml9OmuT0OU+FGhiyTI0k17eGb+ep57p4o9+xs5+die/NdVp3LJmGPIy9ZZKxJ/KnSRONp7oJEXFmzkmfnreK9qN7mRLCaOHsh1Zw1h3LC+Oo9cEkqFLnKUmqLOWxXb+d27G5i5aBO19U2MHFDInX8/in84rUSHVaTLqNBFjoC7s2D9Ll54byMzF25ia80BuudGuOzUY7l67GBOG9xbe+PS5VToIgFFo86763fxx/c389KizazbUUtuJIsLTyzmslNL+LsT+1OQq2PjEh4VukgH9jc08fqqbbz6/hZeW7qVbXsPkJ1lnH1cP275u+P55MkDdTWnJA0Vukgr7s6abfv426ptzFm5jb+t3EZdQxM98rK5YGQxF48awAUj+6vEJSmp0CXj7dxXz+urm8t7zsptbNhVB8CgPgVceWYJF48ayPjhfXWqoSQ9FbpknE2765hfuZPyyh2UV+5k6eY9uENhfjbnHNePmy44jgnHFzG0Xze9sSkpRYUuaa2hKcqKLTW8u24X5ZU7mF+58+AeeLfcCGcM6cO/XHQC544oYkxJL7IjWSEnFjlyKnRJGw1NUVZu2cviDbtZuGEXizbsYemmPdQ3RgEoLsxjXGlfpk4YRtnQvpx0TKEKXNKKCl1SjruzteYAyzbXsGJzDcu31LBiSw3LN9dwIFbePfKyGV3SkxvOKWV0SS9OG9SbwX0LdAhF0poKXZJWfWOU9Ttrqdy2jzWx28ote1m+pYbddQ0HtysuzGPkgEI+O34oYwb14pSSXpT2605WlspbMosKXUJVs7+BDbvq2LirjnXba6ncXsuabfuo3L6Pqp11NEU/mKjVMz+bEwYUcumYYxg5sJATBjTf+urSehFAhS4JtL+hieqaA2ytOcCm3c2lvWFnHRt21VG1s3l5z/7GD31P99wIpUXdGV3Si8tOPZbSft0pLerOsKLu9OmWo0MmIh1Qocth2d/QxK7aBnbW1rNzXz3Vew8cLO2te/Y3f43db1vWAIV52ZT0KaCkdwHjhvXl2N7N90v6FDCoTwHFPfJU2iJHSIWegZqizt79jezZ30DN/kZqWr4eaGBPXeMHhV1bz4599bHybn6str6p3efMy86if888+hfmc3xxD845rh/9C5uXi3vmMbBnPiV9CuiZryssRRIlUKGb2UTgv4EI8Ji739NmvcXWTwZqgRvc/Z04Z81IjU1R6hqaqKtvojZ2q2to/OD+wcdbPxa739DE3taFHbu/7xCl3FphfjZ9uuXSp3suxT3yOGFAIX265dK3e27z491y6N0tl+LCXIoL8+mZn609a5GQdVroZhYBHgIuBqqA+Wb2gru/32qzScCI2O0s4OHY16QSjTpN7jRFnWjsa+tbQ9RpbIrS0OQ0RqM0NjkNTVEao7Gvsccbmvzg/Q8/1rJt7HmizoHGJg40RKlvinKgIcqBxibqG6McaGx7Pxq739TqfvRDbwoGkRMxCnIidMvNpltuhB752RTmZ9O/MJ/C/GwK83NiX7PpefD+B48V5ufQqyCH3Gydny2SaoLsoY8DVrl7BYCZPQ1cDrQu9MuBX7i7A3PNrLeZHePum+Id+K/Lt3LXi+8TdQ4WcfRQJe1ONMrBEg9DbnYWedlZ5GVHYl+zPvRY97xs+nTLIi+neTk30nK/ebvcSISC3CwKcrPplhOhW26Egtzm7yuILXfLzaYgt/l+ji6UEclYQQq9BFjfarmKj+59t7dNCfChQjezacA0gCFDhhxuVgAK83M4cWBPsrKMiBH7akSy7EP3W25ZZkSyIGIfrM+Krcs+uL75lhMxsrOyyI4YOZEssrNiX2OP50SM7FaPH/yeSBY5WbF1ESMnq+V7TIchRKTLBCn09hqp7e5ukG1w9+nAdICysrIj2mU+c2gfzhza50i+VUQkrQX5/bwKGNxqeRCw8Qi2ERGRBApS6POBEWY2zMxygWuAF9ps8wLwOWs2HtidiOPnIiJyaJ0ecnH3RjO7BXiF5tMWZ7j7EjO7Kbb+EWAWzacsrqL5tMUbExdZRETaE+g8dHefRXNpt37skVb3HfhKfKOJiMjh0DluIiJpQoUuIpImVOgiImlChS4ikias+f3MEF7YrBpYG3DzImBbAuPEm/ImXqplVt7EyqS8Q929uL0VoRX64TCzcncvCztHUMqbeKmWWXkTS3mb6ZCLiEiaUKGLiKSJVCn06WEHOEzKm3iplll5E0t5SZFj6CIi0rlU2UMXEZFOpEyhm9kzZrYgdqs0swVhZ+qMmX3VzJab2RIzuy/sPB0xs++Y2YZWP+PJYWcKwsz+zczczIrCztIRM7vLzBbGfrZ/NLNjw87UETO738yWxTL/zsx6h52pI2b2mdi/s6iZJe3ZLmY2MdYJq8zstng/f8oUurtf7e6nuftpwLPAcyFH6pCZXUjzaL4x7n4y8IOQIwXxw5afcewD2ZKamQ2medbturCzBHC/u4+J/f19Efh2yHk68yow2t3HACuA20PO05nFwKeA2WEHOZRW85knAaOAKWY2Kp6vkTKF3sKaZ7pdBTwVdpZO3Azc4+4HANx9a8h50tEPgW/QznSsZOPue1otdifJM7v7H929MbY4l+ahNUnL3Ze6+/Kwc3Ti4Hxmd68HWuYzx03KFTowAdji7ivDDtKJE4AJZvaWmf2fmY0NO1AAt8R+xZ5hZkk958/MLgM2uPt7YWcJyszuNrP1wHUk/x56a18AXgo7RBo41OzluAn0eehdxcxeAwa2s+oOd38+dn8KSbJ33lFemn+2fYDxwFjgN2Y23EM8raiTvA8Dd9G853gX8ADN/5BD00nebwKf6NpEHevs76+73wHcYWa3A7cAd3ZpwDaC/HszszuARuDJrszWnoD9kMwCzV4+GklV6O5+UUfrzSyb5uNkZ3ZNoo51lNfMbgaeixX4PDOL0vz5DdVdla+tzn6+LczsUZqP84bqUHnN7BRgGPBe8xE4BgHvmNk4d9/chRE/JOjPF/g1MJOQCz3Av7fPA5cCHw9zR6TFYfx8k1XCZy+n2iGXi4Bl7l4VdpAAfg/8HYCZnQDkksQfHmRmx7RavILmN5mSkrsvcvf+7l7q7qU0/0M5I8wy74yZjWi1eBmwLKwsQZjZROBW4DJ3rw07T5oIMp/5qCTVHnoA15Akh1sCmAHMMLPFQD3w+WTYy+nAfWZ2Gs2/AlYCXw41Tfq5x8xGAlGaP2X0ppDzdOYnQB7wauy3oLnunrSZzewK4MdAMTDTzBa4+ydDjvUhh5rPHM/X0JWiIiJpItUOuYiIyCGo0EVE0oQKXUQkTajQRUTShApdRCRNqNBFRNKECl1EJE2o0EVE0sT/B5YH8QVmHGIsAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = np.arange(1, 1000)\n",
    "plt.plot(np.log(1/N), np.log(1 + 1/(.35*N)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ad8c660",
   "metadata": {},
   "source": [
    "# Shock to the cost of entry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5526d421",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
